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Chapter  1 


INTRODUCTION 

1.0  Introduction 

In  the  spectral  estimation  of  signals  one  comes  across  a 
problem  when  components  of  the  spectral  response  are  spaced 
closely  together,  and  the  problem  is  further  complicated 
when  the  signal  is  noisy.  There  is  a  point  up  to  which  it 
is  possible  to  resolve  two  pulses  that  are  placed  closely 
together.  Beyond  this  point  it  can  only  be  done  if  some  a 
priori  information  is  given  about  the  signal.  The  point 
at  which  it  becomes  difficult  without  a  priori  knowledge 
can  be  obtained  from  the  sampling  theorem,  and  is  given  by 

T • B=0 . 5  (1.0-1) 

where  the  T  represents  the  duration  of  time  that  the  sig¬ 
nal  is  observed,  and  B  is  the  bandwidth  of  the  signal. 
Equation  (1.0-1)  states  that  the  limiting  time-bandwidth 
product  is  equal  to  0.5.  If  the  product  is  greater  than 
0.5  it  is  possible  to  obtain  a  good  estimate  of  the  spe¬ 
ctrum  ev.en  when  noise  is  present.  On  the  other  hand,  if 
the  product  is  less  than  0.5  it  is  extremely  difficult  to 


resolve  two  peaks  that  are  spaced  closely  together.  The 
only  way  to  be  able  to  resolve  the  two  peaks  is  to  have 
knowledge  about  the  signal.  This  knowledge  can  then  be 
formulated  to  give  one  a  set  of  linear  equations  with 
equality  and  inequality  constraints.  These  equations  along 
with  a  minimizing  procedure  can  be  used  to  obtain  or  im¬ 
prove  the  esimate.  In  this  thesis  we  show  a  method  to 
obtain  an  estimate  of  two  closely  spaced  pulses  in  the 
spectral  domain  when  the  time-bandwidth  product  is  less 
than  0.5*  The  knowledge  that  we  assume  is  that  there  are 
only  two  pulses  in  the  frequency  domain. 

In  signal  processing  there  is  a  class  of  signals  that  one 
comes  across  having  the  property  that  they  only  exist  for 
positive  time.  This  class  of  signals  are  termed  causal. 
If  the  signal  is  also  real,  then  there  is  a  special  re¬ 
lationship  between  the  real  ard  imaginary  parts  of  its 
Fourier  transform.  The  relationship  that  exists  is  that 
the  real  and  imaginary  parts  of  the  Fourier  transform  are 
related  by  the  Hilbert  transform.  From  this  one  need  only 
know  the  real  or  imaginary  part  of  the  spectrum  and  the 
other  part  can  be  easily  calculated.  The  dual  to  this  is 
that  the  spectrum  contains  only  real  components  for  posi¬ 
tive  frequencies.  Then  the  real  and  imaginary  parts  of  the 
time  signal  are  related  through  the  Hilbert  transform. 
This  type  of  time  signal  is  given  the  name  analytic.  One 


of  the  uses  of  the  Hilbert  transform  is  in  the  area  of 
phase  retrieval.  The  problem  can  be  stated  that  given 
only  the  modulus  of  a  signal,  whether  it  be  in  the  time  or 
frequency  domains,  determine  the  phase.  Under  the  con¬ 
ditions  that  the  signals  are  analytic  or  causal,  the 
Hilbert  transform  holds  between  the  modulus  and  phase.  It 
is  a  well  known  fact  that  when  determining  the  phase,  if 
the  signal  has  zeros  that  occur  in  the  upper  half  of  the 
complex  z-plane,  then  the  phase  cannot  be  found.  For  this 
reason  processing  of  the  modulus  has  to  be  performed  before 
or  after  the  estimate  cf  the  phase  has  been  obtained. 


1.1  Spectral  Estimation  of  Two  Pulses 


When  it  is  desired  to  resolve  two  pulses  spaced  closely 
together  with  a  time-bandwidth  product  of  less  than  0.5, 
more  information  has  to  be  known  a  priori.  The  infor¬ 
mation  that  we  are  given  is  that  there  are  two  pulses,  and 
this  is  all  that  is  needed  to  resolve  the  two  pulses.  The 
method  of  solution  is  based  on  knowing  that  two  pulses  in 
the  frequency  domain  corresponds  to  the  sum  of  two  cosines 
in  the  time  domain.  Since  we  are  using  digital  signal 
processing  methods  where  we  are  sampling  the  time  function, 
we  can  represent  the  two  cosines  in  terms  of  their 
Z-transform.  From  the  Z-transform  we  can  obtain  a  set  of 
linear  equations  plus  a  set  of  inequalities  based  on  the 
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coefficients  of  the  Z-transform.  These  equations  are  then 
used  with  mathematical  programming  techniques  to  provide  an 
optimal  transfer  function  to  describe  the  two  pulses  that 
originally  smeared  into  one  when  the  time-bandwidth  prod¬ 
uct  is  less  than  0.5* 


1.2  Modified  Hilbert  Transform 


Let  us  consider  the  evaluation  of  the  Hilbert  transform. 
It  involves  the  convolution  of  a  signal  with  a  kernel 
given  by 


K=l/t . 


(1.2-1) 


If  numerical  methods  are  used  to  evaluate  the  convolution, 
then  difficulties  arise  because  the  kernel  has  a  singu¬ 
larity  which  makes  certain  values  go  to  infinity  during 
its  evaluation.  To  overcome  this  difficulty,  the  fast 
Fourier  transform  algorithm  has  been  used  to  evaluate  the 
convolution.  In  certain  applications  this  leads  to  anoth¬ 
er  problem.  One  may  find  that  the  frequency  spread  of  the 
function  times  the  kernel  exceeds  the  frequency  domain 
defined  by  the  inverse  of  the  sampling  distance,  and  con¬ 
sequently  the  Fourier  transform  is  distorted.  In  this 
way,  the  computer  evaluation  of  the  convolution  integral 
yields  large  errors. 


A  new  form  of  the  Hilbert  transform  is  given  whose  kernel 
does  not  have  the  problems  that  exist  for  the  other  form 
of  the  Hilbert  transform.  This  new  Hilbert  transform  will 
be  applied  to  the  phase  retrieval  problem  to  show  that  simi¬ 
lar  results  can  be  obtained. 


Chapter  2 


ESTIMATION  OF  CLOSELY  SPACED  FREQUENCIES 
BURIED  IN  NOISE  USING  LINEAR  PROGRAMMING 


2.0  Introduction 


In  signal  estimation  one  finds  that  there  exists  a 
time-frequency  duality  in  that  the  signal  may  be  estimated 
either  in  the  time  or  frequency  domain.  If  a  portion  of 
the  signal  is  known  for  a  specific  length  of  time,  the 
method  where  the  signal  is  estimated  outside  the  time  inter¬ 
val  is  referea  to  as  signal  extrapolation.  If  on  the  other 
hand  the  signal  is  observed  over  the  same  length  of  time 
and  from  these  observations  the  power  spectral  density  of 
the  signal  or  simply  the  spectrum  is  to  be  determined,  the 
process  is  called  spectral  estimation.  If  the  signal  is 
known  over  an  infinite  interval,  the  Fourier  transform  of 
its  autocorrelation  yields  the  power  spectral  density.  In 
either  case,  if  the  signal  is  estimated  in  one  domain  the 

estimate  of  the  signal  in  the  other  domain  can  be  obtained 

0 

by  using  the  Fourier  or  Z  transforms. 

There  are  numerous  methods  that  have  been  developed  for 
signal  extrapolation  and  spectral  estimation.  In  extra- 


polating  a  signal  outside  its  observation  interval,  it  is  a 
well  known  fact  that  if  the  signal  is  continuous,  de¬ 
terministic  and  bandlimited  it  is  possible  ^  perform  the 
extrapolation  without  error.  This  is  becav  .  v>  such  a  signal 
is  analytic  and  the  Taylor  series  expansion  can  be  used 
since  in  priciple,  all  the  derivatives  within  the  observed 
interval  of  the  signal  can  be  evaluated.  This  is  called 
analytic  continuation.  However,  in  practice  the  above 
procedure  is  not  feasible  because  if  the  observed  data 
contains  noise,  the  evaluation  of  the  derivatives  would  be 
inaccurate  since  the  derivative  is  a  noise-sensitive  proc¬ 
ess.  Slepian  et  al.  [i]  proposed  an  algorithm  based  on  a 
series  expansion  in  terms  of  basis  functions  called 
prolate  spheroidal  wave  functions.  These  functions  are 
orthogonal  over  the  observation  interval  as  well  as  over 
the  infinite  interval.  The  series  coefficients  can  be 
evaluated  from  the  observations  given  but  the  series  is 
valid  for  all  time.  This  method  is  limited  by  noise  and 
truncation  errors.  In  general  the  numerical  generation  of 
these  functions  is  very  difficult.  A  method  proposed  by 
Papoulis  [2]  reduces  the  mean-square  error  between  the 
estimate  and  the  original  (time-unlimited)  signal  at  suc¬ 
cessive  iterations.  Using  this  error  energy  reduction 
procedure  and  the  properties  of  the  prolate  spheroidal  wave 
functions,  Papoulis  has  shown  that  this  algorithm  converges 
to  the  original  time  unlimited  signal.  The  error  energy 


reduction  procedure  which  Papoulis  uses  was  first  used  by 
Gerchberg  [3]  in  order  to  extend  a  known  segment  of  an 
analytic  spectrum.  •  Gerchberg  viewed  the  limited  spectrum 
as  the  sum  of  the  true  spectrum  and  an  error  spectrum.  In 
the  algorithm  proposed  it  was  not  the  idea  to  construct  the 
spectrum  perfectly  out  to  some  new  limit  but  to  reduce  the 
energy  (squared  function  integrated)  of  the  error  spec¬ 
trum.  The  error  spectrum  is  equal  to  and  opposite  to  the 
true  spectrum  in  the  area  outside  where  the  true  spectrum 
is  known.  The  procedure  was  shown  to  be  very  effective 
against  noisy  data. 

For  discrete- time  signals,  the  analycity  property  vanishes 
due  to  sampling.  Therefore,  the  extrapolated  estimate  need 
not  coincide  with  the  original  signal.  Other  constraints 
besides  the  band-limited  assumption  must  be  imposed  on  the 
estimate  to  achieve  a  unique  solution.  Much  of  the  work  on 
extrapolation  of  discrete-time  signals  is  recent.  Sabri 
and  Steenaart  [4]  suggested  a  discrete  version  of  the 
iterative  algorithm  proposed  by  Papoulis  [2].  The  al¬ 
gorithm  derived  a  solution  by  finding  an  extrapolation 
matrix.  Cadzow  has  proposed  a  different  extrapolation 
matrix  which  does  not  have  the  existance  problems  of  the 
extrapolation  matrix  suggested  in  M  and  has  some  dimen¬ 
sionality  advantages.  Much  of  the  theory  of  extrapolation 
has  been  developed  for  continuous-time  signals,  and  a  solu- 
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tion  for  the  discrete  case  has  been  obtained  by  sampling 
the  continuous-time  solution. 

Estimation  of  the  spectrum  of  discretely  sampled  deter¬ 
ministic  signals  has  usually  been  based  on  procedures  using 
the  fast  Fourier  transform.  These  methods  are  compu¬ 
tationally  efficient  and  produce  reasonable  results  for  a 
large  class  of  applications.  Even  with  these  advantages 
the  fast  Fourier  transform  approach  has  two  limitations. 
The  first  and  most  proninant  limitation  is  the  frequency 
resolution,  i.e.,  the  ability  to  distinguish  the  spectral 
responses  of  two  or  more  signals.  The  frequency  resolution 
is  approximately  given  by  the  reciprocal  of  the  time  inter¬ 
val  over  which  the  sampled  data  is  available.  The  second 
limitation  is  the  windowing  of  the  data  when  processing 
with  the  fast  Fourier  transform.  This  windowing  occurs 
because  in  order  to  use  the  fast  Fourier  transform  the 
signal  is  truncated  at  some*  time  and  assumed  to  be  periodic 
beyond  this  point  [«].  If  the  window  used  is  a  constant 
amplitude  window,  then  its  Fourier  transform  is  the  sam¬ 
pling  function.  When  in  the  spectrum  there  is  an  impulse, 
the  convolution  of  this  impulse  with  the  sampling  function 
produces  a  spectrum  where  the  energy  of  the  impulse  is  no 
longer  localized  at  one  frequency  but  leaks  into  the  side- 
lobes  of  the  sampling  function.  This  leakage  can  have  an 
effect  on  other  spectral  responses  that  are  present.  These 


effects  can  cause  weak  signal  responses  to  be  masked  by 
higher  sidelobes  from  the  stronger  spectral  responses.  By 
selecting  tapered  data  windows  such  as  Hanning,  Hamming  or 
Bartlet  windows  [7]  ,  the  sidelobe  leakage  can  be  reduced 
but  at  the  expense  of  reduced  resolution. 

These  two  limitation  of  the  fast  Fourier  transform  are 
troublesome  when  short  data  records  are  analyzed.  Short 
data  records  are  a  common  occurance  because  many  measured 
processes  are  brief  in  duration.  Some  applications  where 
short  data  records  are  encountered  are  neurophysics  [8,9], 
geophysics  [ lo] ,  speech  communication  [ll,12],  radar  [ 1 3 J • 
sonar  I*].  and  direction  finding  [15].  Due  to  the  inher¬ 
ent  limitations  of  the  fast  Fourier  transform,  many  alterna¬ 
tive  spectral  estimation  procedures  have  been  proposed. 

In  a  classical  paper  by  Blackman  and  Tukey  [l6],  a  prac¬ 
tical  implementation  of  Wiener's  [17]  autocorrelation 
approach  to  power  spectrum  estimation  when  using  sampled 
data  sequences  was  used.  The  method  first  estimates  the 
autocorrelation  lags  from  the  measured  data,  windows  the 
autocorrelation  estimates  in  an  appropriate  manner,  and 
then  Fourier  transforms  the  windowed  lag  estimates  to  ob¬ 
tain  the  power  spectral  estimate.  This  method  was  very 
popular  until  the  introduction  of  the  fast  Fourier  trans¬ 
form  algorithm  by  Cooley  and  Tuckey  £ 1 8 J .  This  renewed 


interest  in  the  periodogram  which  is  obtained  as  the 
squared  magnitude  of  the  output  values  from  a  fast  Fourier 
transform  performed  directly  on  the  data  set.  The  maximum 
entropy  spectral  estimation  method  proposed  by  Burg  [ 1 9 3  is 
based  upon  an  extrapolation  of  a  segment  of  a  known 
autocorrelation  function  for  lags  which  are  not  known.  In 
this  way  the  characteristic  smearing  of  the  estimated  spec¬ 
trum  due  to  the  truncation  of  the  autocorrelation  function 
can  be  removed.  Burg  argued  that  the  extrapolation  should 
be  made  so  that  the  time  series  characterized  by  the  extra¬ 
polated  autocorrelation  function  has  maximum  entropy.  In 
linear  prediction  [20],  the  signal  is  modeled  as  a  linear 
combination  of  past  outputs  and  inputs.  When  only  the 
present  input  is  used,  the  coefficients  of  the  model  are 
solved  for  by  minimizing  the  total  squared  error  with  re¬ 
spect  to  each  of  the  coefficients.  A  very  efficient  al- 
gorihtm  has  been  developed  by  Levinson  [21]  and  Durbin  [22] 
for  solving  the  equation  for  the  coefficients.  A  more 
extensive  overview  is  given  in  papers  by  Jain  et  al.  [23] 
on  extrapolation  algorithms  for  discrete  signals  and  Kay 
et  al.  [24 J  on  spectrum  estimation. 

A  major  concern  of  spectral  and  time  series  estimation  is 
that  of  system  modeling.  Often  if  there  is  more  knowledge 
about  the  process  from  which  the  data  is  taken,  one  is 
able  to  make  a  more  reasonable  assumption  other  than  the 


assumption  that  the  data  is  zero  outside  the  observed  inter¬ 
val.  A  priori  information  or  assumptions  may  permit  selec¬ 
tion  of  an  exact  model  for  the  process,  or  at  least  a  model 
that  is  a  good  approximation  to  the  actual  process.  It  is 
usually  possible  to  obtain  a  better  spectral  estimate 
based  on  the  model  by  determining  the  parameters  of  the 
model  from  the  obsevations.  Ey  using  modeling,  spectral 
estimation  becomes  a  three  step  process.  The  first  step  is 
to  select  a  time  series  model.  The  second  step  is  to 
estimate  the  parameters  of  the  assumed  model  using  the  data 
samples  or  auto  correlation  lags.  The  third  step  is  to  use 
the  estimated  parameters  and  substitute  them  into  the 
theoretical  spectrum  implied  by  the  model.  The  motivation 
for  the  modeling  approach  to  spectral  estimation  is  the 
higher  frequency  resolution  achieved  over  the  traditional 
techniques  such  as  the  periodogram  [24  J .  It  is  easy  to  see 
that  if  one  is  successful  in  developing  a  parametric  model 
for  the  behavior  of  some  signal,  then  the  model  can  be  used 
for  different  applications. 

In  time  series  analysis,  the  continuous-time  signal  s(t) 
is  sampled  to  obtain  a  discrete-time  signal  s(nT),  also 
termed  a  time  series,  where  n  is  an  integer  variable  and  T 
is  the  sampling  interval.  The  sampling  frequency  is  given 
as  f=l/T.  The  expression  s(nT)  can  be  abbreviated  as  s(n) 
by  normalizing  the  discrete  time  scale  by  T. 


r.  r. 


The  most  powerful  model  in  use  today  is  the  ARI«1A  model 
where  a  signal  s(n)  is  considered  to  be  the  output  of  some 
system  with  an  unknown  input  u(n)  such  that  ti.-.  following 
relationship  holds 


s(n)=-  £a(i)s(n-i)+G  £b(l)u(n-l)  »b(0)  =  l 

i»i  i»o 


(2.0-1) 


where  a(i) , l-k-p,b(l) ,  1-1-q ,  and  the  gain  G  are  the  param¬ 
eters  of  the  hypothesized  system.  In  equation  (2.0-1)  the 
output  s(n)  is  a  linear  function  of  past  outputs  and  pre¬ 
sent  and  past  inputs.  The  signal  is  predictable  from  line¬ 
ar  combinations  of  past  outputs  and  inputs,  and  therefore 
it  has  been  termed  linear  prediction. 


Equation  (2.0-1)  can  also  be  expressed  in  the  frequency 
domain  by  taking  the  Z  transform  of  both  sides  of  the  equa¬ 
tion.  If  H(z)  is  the  transfer  function  of  the  system, 
then  we  have  that 


H(z)  =  S(z)/U(z) 


H|?,b(l)z-' 

l  +  ga(i)z"‘ 


(2.0-2) 


where 


S(z)  =  2  s(n)z“" 

B«0 


(2.0-3) 


is  the  Z  transform  of  s(n),  and  U(z)  is  the  Z  transform  of 
u(n).  H(z)  in  (2.0-2)  is  the  general  pole-zero  model.  The 
roots  of  the  numerator  are  the  zeros  and  the  roots  of  the 
denominator  are  the  poles  of  the  model. 


There  are  two  special  cases  of  the  model  which  are  also  of 
interest : 

1)  all-zero  models  a(i)=0,  l£i-p 

2)  all-pole  model:  b(l)=0,  l-l£q. 

In  statistical  literature  the  all- zero  model  is  known  as 
the  moving  average  (MA)  model,  and  the  all-pole  model  is 
known  as  the  autoregressive  (AR)  model.  The  pole- zero 
model  is  then  known  as  the  autoregressive  moving  average 
(ARMA)  model.  In  this  chapter  only  the  pole- zero  model  and 
the  all-pole  model  will  be  considered. 


In  many  applications  one  finds  that  the  spectrum  of  the 
signal  is  composed  of  two  closely  spaced  impulses.  This 
situation  can  appear  in  radar  where  one  of  the  impulses  is 
due  to  a  jamming  signal  and  the  other  is  produced  by  a 
desired  signal.  If  the  energy  of  the  desired  signal  is 
small  compared  to  the  energy  of  the  jamming  signal,  when 

the  time  over  which  the  total  signal  is  observed  is  small 
the  resultant  spectrum  has  the  two  impulses  smeared  togeth- 


er.  From  this  smeared  spectrum  it  is  desirable  to  locate 
the  signal  with  the  smaller  energy.  Kaveh  [25]  and  Cadzow 
[26]  have  proposed  methods  for  estimating  the  parameters  of 
an  ARMA  model  from  the  observed  data  samples.  In  both 
methods  it  is  found  that  the  order  of  the  ARMA  models  has 
to  be  higher  than  for  the  theoretical  case  when  high  resolu¬ 
tion  of  noisy  signals  is  desired.  The  phenomenon  of  the 
smearing  of  the  pulses  into  one  broader  pulse  is  due  to 
the  multiplication  of  the  actual  time  signal  with  a  rectan¬ 
gular  window  to  produce  the  short  time  observations.  In 
the  frequency  domain  the  transform  of  the  short  time  signal 
is  a  convolution  of  the  desired  transform  with  the  trans¬ 
form  of  the  rectangular  window.  If  the  signal  is  concen¬ 
trated  in  a  narrow  bandwidth,  this  convolution  operation 
will  spread  the  energy  of  the  process  into  adjacent  fre¬ 
quency  regions.  The  convolution  of  the  window  transform 
with  that  of  the  actual  signal  transform  means  that  the 
most  narrow  spectral  response  of  the  resultant  transform 
is  limited  to  that  of  the  main-lobe  width  of  the  window 
transform,  independent  of  the  data.  For  the  rectangular 
window,  the  main-lobe  width  is  approximately  the  inverse  of 
the  observation  time.  When  additive  noise  is  also  pre¬ 
sent,  the  spectral  response  whose  energy  level  is  smaller 
than  the  noise  energy  level  can  be  hidden  by  the  noise 
making  it  difficult  to  determine  that  particular  spectral 
response.  Gerchberg  [3] »  and  Mammone  and  Eichmann  [27»28] 


have  devised  methods  for  estimating  the  spectrum  when  there 
is  a  lack  of  spectral  resolution  due  to  windowing. 

In  this  chapter  a  new  method  for  estimating  the  coef¬ 
ficients  of  an  ARMA  process  is  introduced.  The  ARMA  proc¬ 
ess  considered  is  one  whose  spectrum  contains  two  impulses 
closely  spaced  together  in  frequency  and  the  power  of  one 
impulse  is  larger  than  the  other  impulse.  This  is  the  case 
that  would  be  found  in  radar.  The  data  is  observed  for  a 
very  short  time,  which  along  with  additive  noise  makes  it 
difficult  to  distinguish  the  two  impulses  from  one  another. 
The  method  of  solution  assumes  two  things.  First,  it  is 
assumed  that  the  approximate  region  in  the  frequency  do¬ 
main  where  the  two  impulses  occur  is  known.  The  second 
assumption  is  that  the  signal  is  known  to  be  the  sum  of  two 
sinusoids  whose  spectrum  is  given  by  two  impulses.  This 
assumption  can  be  made  by  looking  at  the  resultant  spectrum 
of  the  time  windov/ed  signal.  From  the  first  assumption  it 
is  known  where  approximately  the  pulses  are  to  be.  In  that 
region  the  spectrum  is  very  broad  which  could  be  only 
possible  if  the  sampling  function  is  convolved  with  two 
impulses  and  not  with  one  impulse.  The  method  of  solution 
uses  linear  programming  to  solve  for  the  coefficients  of 
the  ARMA  model.  Linear  programming  has  been  used  for  spec¬ 
trum  estimation  by  Mammone  and  Eichmann  [27.28]  as  well  as 
Levy  et  al.  [29].  In  this  method  ,  linear  programming  is 


used  to  minimize  the  1,  norm  of  the  absolute  values  of  the 


error.  The  error  is  the  difference  between  the  actual 
spectrum  and  the  one  obtained  from  the  ARMA  model  used  in 
the  region  where  the  response  is  known  to  be.  This  error 
is  sometimes  termed  the  residual  error.  The  new  method 
which  estimates  the  spectrum  and  from  which  the  time  signal 
estimate  can  be  obtained  is  demonstrated  by  computer  simu¬ 
lations  . 


2.1  Parameter  Estimation 


The  method  by  which  the  parameters  of  the  pole-zero  model 
will  be  estimated  is  a  two  step  process,  first,  to  estimate 
the  denominator  coefficients  and  second,  to  estimate  the 
coefficients  of  the  numerator  independently.  V/ hat  this 
method  is  doing  is  to  estimate  the  coefficients  of  the 
transfer  function  H(z).  The  Z  transform  of  the  signal  S(z) 
can  be  thought  of  as  the  output  of  the  filter  H(z),  with 
an  input  function  U(z).  If  the  input  u(n)  is  taken  as  a 
unit  impulse,  then  one  has  that  the  transfer  function  H(z), 
and  the  signal  S(z)  are  identical,  or 


H( z ) =S( z ) . 


(2.1-1) 


Therefore,  by  estimating  the  parameters  of  H(z)  one  is  also 
estimating  the  parameters  of  S(z). 


>>>>>-• 


The  first  step  in  the  estimation  of  the  parameters  will  be 
to  estimate  the  coefficients  of  the  denominator,  a(i). 
Equation  (2.0-1)  can  be  rewritten  in  the  form 


s(n)=-  £  a(i)s(n-i)+B(n) 
i=i  . 


(2.1-2) 


where 

q 

B(n)=G  X  b(l)u(n-l),  b(0)=l. 

1=0 

If  in  (2.1-2)  one  takes  the  summation  to  be  an  approxima¬ 
tion  of  the  signal  s(n),  one  has  that  the  approximation 
s'(n)  is 

s' (n)=-  £  a(i)s(n-i) .  (2.1-3) 

i=i 

Then  the  error  between  the  actual  value  s(n)  and  the  esti¬ 
mated  value  s'(n)  is  given  by 

B(n)=s(n)-s'(n)=s(n)+  £  a(i)s(n-i).  (2.1-4) 

i=i 

The  value  B(n)  is  then  taken  to  be  a  form  of  the  residual. 


This  method  estimates  the  parameters  from  data  in  the  fre¬ 
quency  domain,  so  that  (2.1-2)  is  converted  into  the  fre¬ 
quency  domain  by  taking  the  discrete  Fourier  transform 
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(DFT)  of  both  sides  of  the  equation.  Since  the  fast 
Fourier  transform  (FFT)  will  be  used  in  all  the  simula¬ 
tions,  the  DFT  is  used  instead  of  th^  7  transfo^~  Tvi^ 
is  possible  because  a  periodic  sequence  ..as  a  DFT  which  can 
be  interpreted  as  samples  on  the  unit  circle,  equally 
spaced  in  angle,  of  the  Z  transform  of  one  period.  If  a 
signal  is  periodic  with  period  N,  it  is  possible  to  repre¬ 
sent  this  signal  in  terms  of  a  Fourier  series  consisting  of 
a  sum  of  sines  and  cosines  or  equivalently  complex  expo¬ 
nential  sequences  with  frequencies  that  are  integer  multi¬ 
ples  of  the  fundamental  frequency  2tt/N  associated  with  the 
periodic  sequence  [6] .  The  DFT  analysis  and  synthesis  pair 
are  expressed  as 

X  ( k )  =  £  x ( n ) exp ( -  j 2-irkn/N ) 

n=0 

x(n)  =  l  ^  X(k)exp(j2irkn/N) 

—  k  =  0 

N 

where  both  x(n)  and  X(k)  are  periodic.  This  gives  us 
^  s(n)exp(- j2imk/N)  = 

jT*°  N-l 

-  Y  exp(- j2irik/N )a(i)  Y  s(n)exp(- j2irnk/N) 

2*1  n»0 

+B(exp((j2imk/N)  (2.1-5) 


Equation  (2.1-5)  represents  the  relationship  between  the 
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DFT  of  the  signal  s(n),  the  DFT  of  the  convolution  sum 
(where  p<N,  so  that  the  rest  of  the  a(i)'s  are  set  to 
zero),  and  the  DFT  of  the  residual.  This  equation  can  be 
written  in  a  simpler  form  as 


S(k)=-  [  £a(i)exp(- j2irik/N)]  S(k)+B(k)  (2.1-6) 

where  k  ranges  between  the  values  of  0  and  N-l.  The  var¬ 
iables  S(k)  and  B(k)  are  complex  numbers  so  that  they  may 
be  written  in  the  form 


S(k)=S' (k)+jS* ’ (k) 

B(k)=B'(k)+jB' *(k).  (2.1-7) 


P 

P 


Ey  making  use  of  Euler's  formula  for  the  exponent,  (2.1-6) 
is  written  as 


S' (k)+jS* ' (k)= 

-  f  a(i)  (cos(2trik/N)- jsin(27rik/N) )  (S'  (k)  +  jS"  (k)) 

i*i 

+  (B'(k)  +  jB"(k)) 

=-  H  a(i)  (S' (k)cos(27rik/N)+S' '  (k)sin(2irik/N)  )+B*  (k) 

p 

-  j  (  a(i)  (S' '  (k)cos(2Trik/N)-S'  (k)sin(27rik/N)  )-B'  '  (k) ) ) 

(2.1-8) 


From  (2.1-8)  there  are  two  sets  of  equations,  one  for  the 
real  part  and  the  other  for  the  imaginary  part  given  by 


w 


■  -  •  .A  - 


■■ 


s*  (k)=-  |a(i)(S'  (k)cos(27rik/N)  +  Sl  1  (k)sin(2irik/N)  )+B'  (k) 

(2.1-9) 

S'  ’  (k)=-  f  a(i)  (S’  ’  (k)cos(27rik/N)-S'  (k)sin(27rik/N)  )+B'  1  (k) 

i=i 

(2.1-10) 

Eecause  linear  programming  will  be  used  to  solve  for  the 
coefficients  a(i),  we  can  also  specify  the  range  of  the 
a(i)'s  to  help  in  the  solution.  The  range  of  the  coef¬ 
ficients  a(i)  can  be  specified  because  it  is  known  that 
the  signal  s(n)  is  made  up  of  the  sum  of  twro  cosines  of 
different  but  closely  spaced  frequencies.  The  denominator 
of  the  Z  transform  of  a  cosine  wave  can  be  derived  as  fol¬ 
lows.  From  the  definition  of  the  Z  transform  and  the  ex¬ 
ponential  representation  of  the  cosine  one  has 

CO  CO 

cos(wkT)  *—  Z 0.5exp( jwkT)z  * +  Z 0.5exp( jwkT) z  k  (2.1-11) 

k*0  k  »  O 

In  each  of  the  summations  the  exponential  can  be  multi¬ 
plied  by  z  and  can  be  written  as 

exp ( jwkT  )z"k  =  (z'1  exp (  jwkT)  )k 

exp  (- jwkT)  z"k  =  (z-1  exp  (- jwkT  ))k  (2.1.12) 

By  substituting  each  term  into  the  appropriate  summation 
and  using  the  infinite  summation  identity  that 


Z a  = 

M»0 


(2.1-13) 


Equation  (2.1-11)  is  written  as 


cos(wkT)  ••  0.5 


[ - i -  * - i - 1 

ll-z-1  exp(jwT)  1-z"1  exp(-jwT). 


=  0.5  1-z"1  exp{- jwT)  +  l-z“'  exp(jwt)  N(z) 

( l-z-'  exp(  jwT) )  ( l-z"'  exp(-jwT))  D(z) 


(2.1-14) 


By  multiplying  the  denominator  one  obtains 


D(z)=  l-(exp(  jwT)+exp(-jwT)  )z“'  +z_1 


(2.1-15) 


Equation  (2.1-15)  can  also  be  reduced  to  the  form 


D(z)  =  l-2coswTz"'+z‘ 


(2.1-11) 


where  wT  is  the  normalized  radian  frequency  of  the  cosine 


wave,  and  T  is  the  sampling  interval.  Knowing  that  s(n) 


has  two  cosine  waves,  the  Z  transform  S(z)  has  to  have  the 


product  of  the  denominator  of  each  of  the  cosine  waves  to 


give 


D'  (z)  =  l-2(cosw,  T+cosw2T)  z~'+2  ( l+2cosw,  TcoswaT)z** 


-2(cosw,  T+cosw2T)  z~s+z~A 


(2.1-17) 


where  w  and  w  are  the  two  radian  frequencies  of  the  co¬ 
sines.  Equation  (2.1-17)  can  be  written  as 

D*  (z)  =  l+a(l)z"'+a(2)  z‘2+a  (  3 )  f®+a  ( 4 )  z"*.  (2.1-16) 

where  a( 1 )=-2(coswf T+cosw2T) 

a(2)=2(  l+2coswt  Tcosw2T) 
a(3)=-2(cosw,  T+cosw2T) 
a(4)=l . 

Each  of  the  coefficients  in  (2.1- IS)  contains  cosines.  The 
limit  of  a  cosine  is  -1  and  +1.  These  upper  and  lower 
limits  of  the  cosine  can  be  used  in  the  expressions  for  the 
coefficients  to  give  the  limits  on  these  coefficients. 
This  gives  the  limits  as 


-4-a( 1 )-k 
-2-a( 2 )-6 

a(3)=a(l) 

a(4)  =  l .  (2.1-19) 

A  system  is  said  to  be  stable  if  all  the  poles  of  its  trans¬ 
fer  function  are  contained  inside  the  unit  circle  of  the 
z-plane.  In  the  special  case  of  sinusoids  one  finds  the 
poles  to  be  located  on  the  unit  circle  so  that  the  system 


Mv.'iiviv'.v(vy^'iVAV 


can  be  termed  marginally  stable.  By  marginally  stable  it 
is  meant  that  the  response  is  oscillating  but  the  poles  are 
not  outside  the  unit  circle.  Since  the  sinusoid  is  an 
oscillating  function,  its  poles  have  to  be  located  on  the 
unit  circle.  Therefore,  if  a  sinusoidal  response  is  termed 
stable,  it  is  only  because  it  has  no  poles  that  are  lo¬ 
cated  outside  the  unit  circle.  If  the  solution  gives  values 
according  to  (2.1-19)  then  one  is  assured  a  stable  filter 
because  one  obtains  these  limits  through  the  Z  transform 
of  the  cosine  wave  which  can  be  called  a  stable  function. 

Equation  (2.1-6)  can  be  written  as 

(1+  X  a(i)exp(-j27rik/N))S(k)=B(k)  (2.1-20) 

*»• 

By  evaluating  the  z  in  (2.1-16)  on  the  unit  circle  one  has 
D ' ( exp ( j 2irik/li ) )  =  1  +  £  a(i)exp(- j27rik/N)  (2.1-21) 

•*l 

which  is  seen  to  be  identical  to  the  term  in  the  brackets 
of  (2.1-20)  with  p=4.  Therefore,  we  know  that  there  are 
only  four  unknown  coefficients  in  the  denominator  of  the 
filter  H(z).  Based  on  the  assumption  that  the  signal  con¬ 
tains  tv/o  sinusoids,  it  is  possible  to  determine  the  exact 
number  of  coefficients  for  the  denominator  of  the  filter 
H(z).  This  is  possible  because  the  method  of  solution 


will  be  able  to  give  results  which  are  close  to  the  actual 
values  even  when  noise  is  added.  In  methods  where  noise 
plays  a  large  role  in  the  estimation  of  coefficients,  the 
order  of  the  models  which  determine  the  number  of  coef¬ 
ficients  used  has  to  be  larger  than  the  theoretical  number. 
In  many  such  cases  the  choice  of  the  model  order  is  based 
on  trying  different  orders  until  the  best  results  are 
obtained  [27]. 


If  one  was  only  interested  in  obtaining  the  frequencies  of 
the  signal,  by  factoring  the  denominator  polynomial  one 
could  obtain  the  radian  frequencies  of  the  signal  from  the 


argument  of  the  complex  roots.  Since  it  is  also  of 
interest  to  obtain  the  time  signal,  it  becomes  necessary 
to  go  on  to  the  second  step  of  the  estimation  process  where 
one  has  to  first  do  some  rearranging  of  (2.0-2)  in  order  to 
estimate  the  coefficients  of  the  numerator  in  a  manner 
similar  to  the  one  just  described.  It  is  not  possible  to 
use  B(k)  that  comes  from  the  solution  of  the  first  part  for 
the  coefficients  of  the  denominator,  because  if  the  signal 
has  noise  added  onto  it,  by  minimizing  this  residual,  the 
noise  becomes  part  of  this  residual  term.  Therefore,  the 


first  thing  one  has  to  do,  is  to  take  (2.0-2)  and  invert 


it,  to  give 


G 


(2.1-22) 


=  1+  ^  a  (  i )  z~‘ 

H(z)  1+  p(l)r' 

v;here  for  simplicity  we  will  replace  the  H(z)  tern  by  S(z) 
and  the  quotion  G/S(z)  by  S*(z).  This  is  done  because  the 
sane  method  is  used  to  estimate  the  coefficients  b(l). 
The  gain  G  cannot  be  found  using  this  method  so  that  the 
signal  estimate  will  always  be  normalized  to  one.  From 
(2.1-22)  we  obtain  an  expression  similar  to  (2.1-2)  with  a 
residual  A(n)  given  by 


s*(n)=-  Z  b(l)s*(n-l)+A(n) 

!=l 

where  s*(n)  is  the  inverse  Z  transform 
through  the  sane  analysis  as  before  we 
lov/ing  two  equations  similar  to  (2.1-9)  and 


( 

of  S*  (  z ) 
obtain  t 
(2.1-10) 


he 


I 


1-23) 


Going 

fol- 


S*’(k)=-  Z  b ( 1) ( S* '  (k)cos( 27rlk/H )  +  S* '  ’  (k)sin(27rlk/N)  )+A*  (k) 

l=*f 

(2. 1-24) 

S*'  '  (k)=-  f)  b(l)  (S*'  '  (k)cos(27Tlk/N) 

+S*’  (k)sin(2  7r  lk/N)  )+A*  •  (k) 

(2.1-25) 


where  S*'(k),A'(k)  are  the  real  parts  and  S* 
are  the  imaginary  parts  of  S*(k)  and  A(k). 


(k).A’ ' (k) 


a 


In  order  to  obtain  the  range  of  values  of  the  b(l)'s,  the 
numerator  N(s)  of  the  Z  transform  of  the  sum  of  tv.’o  cosines 
is  given  by 

N(s)  =  (A,  +AZ)-(A,  (cosv;,  T+2cosw2  T  )+A2  (  cosv.2T  +  2cosw,  ?) )  z~x 
+  (A,  +A2)  ( l  +  2cosw,  Tcosw2T)z"2 

-(A,  cosv;,  T+A2cosv/2T)  s-3  (2.1-26) 

where  A,  and  A2  are  the  amplitudes  of  the  cosines  and  w,  T 
and  w£T  are  as  before.  Again  we  have  that  the  values  of 
the  cosines  in  (2.1-26)  are  between  -1  and  1  so  that  we  can 
obtain  the  range  of  the  b(l)'s.  YJhen  substituting  in  the 
limiting  values  of  the  cosines  into  (2.1-26)  it  is  noticed 
that  each  coefficient  has  the  sun  of  (A, +AS).  This  sun  is 
the  gain  factor  G  which  cannot  be  found  using  this  method 
so  that  it  is  factored  out  of  the  numerator.  Equation 
(2.1-26)  can  then  be  represented  as 

N •  ( z ) *  1+b '  ( 1 ) i*  +b *  ( 2 ) z“£+b *  ( 3 )  z"3  (2.1-27) 

where 

N* (z)*N(c)/(A,  +Aa). 

The  modified  coefficients  of  the  numerator  b'(l)  have  the 
following  ranges 


-3-b* ( l )-3 
-3-b' (2)-3 


-l*b*(3)-l  v  2 . 1-23) 

Again  we  are  assured  a  stable  filter  if  the  coefficients 
are  within  these  ranges.  It  can  be  shown  as  before  that 
(2.1-24)  and  (2.1-25)  can  be  used  to  find  the  numerator 
coefficients  of  H(z)  with  q=3* 

Once  the  coefficients  of  both  the  numerator  and  de¬ 
nominator  have  been  calculated,  one  can  obtain  the  spectrum 
from  H(z)  or  the  normalised  estimate  of  the  time  series  by 
taking  the  inverse  Z  transform  of  H(z).  The  inverse  Z 
transform  is  chosen  for  determining  the  estimated  time 
signal  because  it  is  possible  to  obtain  a  closed  form  solu¬ 
tion  of  the  estimate.  From  this  closed  form  the  ratio  of 
the  amplitudes  and  the  phases  can  be  obtained  to  compare 
with  the  exact  solution.  So  far  there  has  been  no  mention 
of  exactly  how  the  coefficients  are  to  be  solved  for.  In 
the  next  section  a  discussion  of  how  the  problem  is  formu¬ 
lated  to  linear  programming  which  is  the  method  of  solution 
used  to  obtain  the  unknown  values  of  the  coefficients  is 
discussed . 


A  linear  optimization  problem  can  be  stated  in  the  form 
minimise  £  (c(i)-Dx) 

I 

subject  to  Px-q 

Gx=h  (2.2-1) 

v;here  the  term  in  the  summation  is  called  the  objective, 
and  the  two  sets  of  linear  equations  are  the  conditions  of 
the  problem.  An  optimum  solution  is  obtained  when  the 
conditions  of  the  problem  and  the  given  objective  are 
satisfied  simultaneously.  A  minimum  feasible  solution 
satisfies  (2.2-1)  and  the  condition  that  all  the  x's  of  the 
solution  are  nonnegative.  The  minimum  feasible  solution  is 
obtained  using  linear  programming  methods  which  were  first 
introduced  by  Dantzig  [30]* 

In  many  applications  of  the  linear  optimization  problem, 
the  solution  vector  must  have  both  positive  and  negative 
values.  To  overcome  this  problem  the  solution  vector  is 
replaced  by  the  difference  of  two  positive  valued  vectors 
given  by  x-y.  Looking  at  (2.2-1)  a  special  form  of  the 

optimization  is  obtained  if  the  variables  c(i)  and  D  of  the 
objective  are  replaced  by  the  variables  h  and  G  re¬ 
spectively.  It  is  seen  that  by  doing  this  substitution 


the  objective  is  the  summation  of  the  difference  between  a 
value  obtained  by  multiplying  the  solution  vector  x  with 
the  matrix  G,  and  the  known  quantity  h.  This  summation 
which  can  be  termed  the  residual  is  zero  only  when  the 
conditions  are  satisfied  exactly.  In  the  cases  where  they 
are  not,  the  residual  will  be  a  positive  or  negative 
number.  Therefore,  one  can  write  the  objective  as 

minimize  V  (h-Gx+Gy)  (2.2-2) 

or 

minimize  £  (u-v)  (2.2-3) 

where  u  and  v  are  positive  valued  vectors  to  be  consistent 
with  the  linear  programming  formulation.  From  (2.2-2) 
and  (2.2-3)  one  has 

Gx-Gy+u-v=h  (2.2-4) 

where  x-0 
y-o 
u-0 
v^O. 

Using  (2.2-3)  and  (2.2-4),  equation  (2.2-1)  is  rewritten  as 


minimize 


I  u-v 


subject  to  Px-Py  -q 

Gx-Gy+u-v=h  (2.2-5) 

In  this  form  the  linear  optimization  prob'1''-  xs  minimizing 
the  1,  norm  of  the  residual  (u-\  ■ 

The  vai'iables  x,y.  .w  v  of  the  linear  optimization  prob¬ 
lem  can  be  solved  for  using  the  simplex  method  [31].  In 
the  simplex  method  once  a  feasible  solution  has  been  deter¬ 
mined,  a  minimum  feasible  solution  is  obtained  in  a  finite 
number  of  steps.  These  steps,  consist  of  finding  a  new 
feasible  solution  whose  corresponding  value  of  the  objec¬ 
tive  function  is  less  than  the  value  of  the  objective  fun¬ 
ction  in  the  preceding  case.  This  process  continues  until 
a  minimum  solution  has  been  reached.  One  is  never  guaran¬ 
teed  that  a  solution  exists.  If  no  solution  exists  it  is 
either  because  no  solution  exists  in  terms  of  nonnegative 
values  of  the  variables  can  be  found  or  a  ncnnegative  solu¬ 
tion  yields  an  infinite  value  to  the  objective  function. 

It  is  found  that  for  the  problem  of  section  2.1,  it  is 
possible  to  obtain  a  minimum  feasible  solution  in  deter¬ 
mining  the  coefficients  of  the  ARMA  model.  The  next  step 

is  to  show  how  one  goes  about  relating  the  problem  of 
section  2.1  to  the  linear  optimization  problem. 


In  (2.2-5)  it  is  seen  that  the  variable  q  is  part  of  the 
inequality  condition.  When  the  equation  for  the  deter¬ 
mination  of  the  denominator  coefficients  was  being  de¬ 
rived,  it  was  seen  from  (2.1-9)  that  the  coefficients  have 
limiting  values.  Two  of  these  conditions  can  be  placed 
into  the  inequality  condition  of  (2.2-5).  As  an  example 
the  inequality 

.4  *  a( 1 )  -4  (2.2-6) 

can  be  written  as 

a(l)  ^  k 

-a(l)  *  4  (2.2-7) 

The  first  inequality  of  (2.2-7)  represents  the  upper  limit 
and  the  second  inequality  represents  the  lower  limit  of 
(2.2-6).  The  last  two  conditions  of  (2.1-19)  can  be  placed 
into  the  equality  condition  cf  (2.2-5).  Equations  (2.1-9) 
and  (2.1-10)  can  also  be  placed  into  the  equality  condi¬ 
tions  of  (2.2-5).  The  a(i)'s  of  the  equations  are  the 
unknowns  being  solved  for,  the  S'(k)  and  S*'(k)  on  the  left 
side  of  the  equals  sign  make  up  the  vector  h,  the  B’(k) 
and  B’’(k)  are  the  residuals  which  are  to  be  minimised,  and 
the  matrix  G  is  made  up  of  the  terms  involving  the  sines 
and  cosines.  By  making  these  substitutions  into  the  lin- 
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ear  programming  format,  one  is  able  to  solve  for  the  un¬ 
known  coefficients.  An  identical  analysis  can  be  made  for 
determining  the  coefficients  of  the  numerator  from  (2.1-24), 
(2.1-25)  and  (2.1-28). 

2.3  Computer  Simulation  and  Results 

The  method  of  solution  was  simulated  using  the  Simplex 
Linear  Programming  subroutines  contained  in  the  Inter¬ 
national  Mathematical  and  Statistical  Library  (IMSL) 
Fortran  callable  subroutine  package.  CUNY's  IBM  computer 
system  was  used.  Tc  show  how  effective  this  new  method  is, 
we  chose  a  signal  whose  FT  contains  two  frequency  pulses  in 
the  frequency  domain  which  could  be  placed  closely  to¬ 
gether.  To  satisfy  this  condition  we  chose 

s(n)=  y/2.  cos(2  rr  ( 110)n/N)+ *^20  cos( 2 it  ( 1 14 )n/N)  +  w(n) 

(2.3-1) 

where  C-n-511,  N=512  and  w(n)  is  Gaussian  noise  with  var¬ 
iance  c r2.  The  signal-to-noise  ratio  (SNR)  of  the  first 
cosine  is  given  by  10iog(l/<r2)  and  of  the  second 
101og(  10/er2 ) .  The  bandwidth  is  given  by 

Al  =  (l  14-110) /K =0.0  07  8125 
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which  will  be  used  in  the  determination  of  the 
time-bandwidth  product. 

Fig.  2.3.1  shows  the  signal  s(n)  when  there  is  no  noise  and 
all  512  points  are  used.  We  can  see  that  s(n)  contains  an 
envelope  of  a  sinusoid  with  six  periods  in  the  512  points 
shown.  Fig.  2.3*2  shows  the  log  of  the  magnitude  of  the 
Fourier  transform  of  s(n),  where  we  see  tv/o  pulses  located 
at  the  frequencies  of  the  signal.  The  rest  of  the  mag¬ 
nitudes  are  so  close  to  zero  that  their  logs  are  very 
negative  and  cannot  be  shown  in  the  figure.  In  Fig.  2.3.3 
we  have  added  noise  with  a  variance  cr  =1.0  and  we  see  that 
it  is  impossible  to  distinguish  how  many  periods  the  en¬ 
velope  cf  the  signal  contains.  The  spectrum  of  this  noisy 
signal  is  shown  in  Fig.  2.3*^  where  one  of  the  things  that 
one  notices  right  away  is  that  the  spectrum  is  now  centered 
about  the  0.0  db  level.  This  corresponds  to  the  spectrum 
of  the  noise  being  added  to  the  spectrum  of  the  signal. 
The  reason  why  the  noise  did  not  affect  the  pulses  is  be¬ 
cause  the  noise  is  not  large  enough  to  bury  any  of  the 
pulses  when  all  of  the  points  are  used  in  determining  the 
spectrum.  Noise  plays  a  large  role  in  the  determination  of 
the  lcacation  of  these  pulses  when  not  all  the  points  are 
used  and  the  tv/o  pulses  begin  to  smear  together. 


OJi> 


and  zeros  of  s(n),  Fig.  2.3*5  shows  the  first  quadrant  of 
the  z-plane.  l.Te  chose  to  show  only  the  first  quadrant 
because  all  the  poles  and  zeros  of  our  s(n)  can  be  shown 
here  with  the  understanding  that  for  complex  poles  and 
zeros,  there  are  mirror  images  of  these  poles  and  zeros  in 
the  fourth  quadrant.  Fig.  2.3.5  shows  two  poles  located  on 
the  unit  circle  which  represent  the  locations  of  the  two 
pulses  in  the  spectrum.  The  zero  in  between  the  pcles  on 
the  unit  circle  is  the  dip  that  occurs  between  the  two 
pulses  in  the  spectrum  of  Fig.  2.3*2.  The  zeros  on  the 
real  axis  along  with  the  complex  poles  and  zeros  are  used 
to  calculate  the  amplitudes  and  phases  of  the  signal  s(n). 


In  order  to  have  a  uniform  method  of  measure  when  truncat¬ 
ing  the  signal  s(n),  we  define  the  time  bandwidth  product 
(TBW)  as  the  number  of  samples  used  multiplied  by  the  nor¬ 
malized  bandwidth  of  the  signal  which  was  defined  earlier 
for  our  case.  This  definition  can  be  explained  by  looking 
at  the  spectrum  of  a  DFT  having  length  N .  The  DFT  is  a 
two-sided  transform  so  that  the  spectrum  has  both  positive 
and  negative  frequencies.  The  maximum,  frequency  of  the 
spectrum  is  N/2  which  in  terms  of  a  normalized  frequency 
can  be  divided  by  N  to  give  1/2.  This  frequency  is  known 
as  the  Nyquist  frequency  because  it  gives  the  separation 
between  time  samples  as  l/N.  The  Nyquist  frequency  is  also 
the  bandwidth  of  the  signal  since  it  is  assumed  that  the 


spectrum  has  a  periodic  extension  beyond  this  frequency. 
Therefore,  if  all  N  time  samples  are  used,  the 

tine-bandwidth  product  is  given  by 
T3V;=  N-  (1/H)  •  (N/2)  =  N/2 

from  which  the  normalized  TEW  is  given  by  1/2.  For  the 
DFT,  the  separation  of  tine  samples  is  given  by  l/N.  If 
the  bandv/idth  of  the  signal  is  taken  to  be  less  than  N/2, 
it  is  seen  that  the  normalized  bandwidth  will  be  less  even 
when  all  of  the  N  samples  are  used.  A  normalized  TBW  of 
less  than  0.5  corresponds  to  choosing  the  separation  be¬ 
tween  time  samples  too  far  apart  so  that  the  signal  is 
undersampled.  It  can  be  seen  that  a  TBW=0.5  is  the  limiting 
value  for  spectral  resolution.  If  TEW  becomes  smaller,  then 
any  two  peaks  spaced  closely  together  will  become  un- 
distinguishable  because  they  smear  into  one  broader  peak. 
Even  when  TBW  is  greater  than  0.5*  there  may  be  difficulty 
in  distinguishing  the  two  peaks  if  the  noise  power  is  high 
enough  to  hide  the  weaker  pulse. 

For  the  first  case  to  test  our  method,  we  chose  to  use  only 
92  out  of  the  512  samples.  This  gives  a  TBW  of  O.719.  In 
Fig.  2.3.6  we  show  the  truncated  time  signal  s(n)  with  noise 
variance  <r2=1.0.  The  spectrum  obtained  using  the  FT  of 
this  signal  is  shown  in  Fig.  2.3.7.  In  the  spectrum  of 
Fig.  2.3.7  we  see  that  the  two  pulses  are  smearing  into  one 
with  a  very  shallow  dip  between  the  two  peaks  very  close 
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in  height  to  the  weaker  pulse.  This  shallow  dip  is  due  to 
both  the  TBW  and  the  noise.  The  combination  of  the  two 
makes  for  the  difficulty  in  judging  that  there  should  be 
two  peaks.  After  using  our  method,  in  calculating  the 
coefficients  of  the  numerator  and  denominator  of  the  para¬ 
metric  model,  we  found  the  pulses  of  the  spectrum  to  be  at 
normalized  frequencies  f,  =0.2148  and  f^ =0.2285  as  compared 
to  the  actual  values  f,  =0.2149  and  f2  =0.222?  which  shows 
that  we  are  very  close  to  the  original  values  with  our 
estimates.  The  pulses  at  these  frequencies  are  shown  in 
Fig.  2.3*8  where  the  dip  that  is  located  between  the  two 
pulses  is  at  f =0.2231  as  compared  to  where  it  should  be  at 
f=0.2l6?.  The  location  of  the  dip  represents  the  complex 
zero  located  on  the  unit  circle  of  the  z-plane. 


We  observe  that  the  first  pulse  in  Fig.  2.3.8  is  larger 
than  the  second  one,  and  if  we  look  at  Fig.  2.3.2  we  see 
the  opposite.  If  we  take  a  look  at  the  location  of  the 
poles  and  zeros  from  the  coefficients  calculated  we  see 
that  they  are  located  as  follows 

poles:  0.219+j0. 975  zeros:  0.l66+j0.977 

0.134+j0.991  0.203 

0.0 

which  are  plotted  in  the  z-plane  in  Fig.  2.3*9*  The  actual 
locations  are  shown  in  Fig.  2.3.5  and  are  given  by 
poles:  0. 219+j0.975  zeros:  0.20?+j0,978 


VJV 


0. 170+J0.9S5 


0.183 

0.0 


By  comparing  the  location  of  the  poles  and  zeros  we  see 
that  there  is  an  error  in  the  location  of  the  estimated 
ones.  These  errors  are  not  that  great  but  they  have  caused 
the  first  peak  to  be  higher  than  the  second  peak.  By  tak¬ 
ing  the  inverse  2  transform  of  the  parametric  model  using 
the  coefficients  calculated,  we  obtained  the  normalized 
signal  estimate  s(n)  given  by 

s(n)=0.6l2cos(1.3^9n+0. 1479) +0.402cos(1.435n-0. l?4l) 

(2.3-2) 

which  if  we  compare  to  (2.3-1)  we  see  that  the  first 
cosine  has  a  larger  amplitude  and  both  cosines  have  a  phase 
which  is  not  in  (2.3-1).  The  estimated  signal  is  shown  in 
Fig.  2.3*10  where  the  envelope  of  the  signal  is  a  sinusoid 
having  seven  periods  in  the  512  samples  compared  to  six  for 
the  original  signal.  Again  this  is  a  consequence  of  the 
error  of  the  location  of  poles  and  zeros. 

For  the  next  case  investigated  we  chose  the  same  TBW  prod¬ 
uct  as  before  but  increased  the  noise  power  too-2  =5*0.  The 


individual  SNR'S  are  given  by  -7  db  for  the  smaller 
amplitude  cosine  and  3  db  for  the  larger  amplitude  cosine. 
Fig.  2.3*11  shows  the  time  truncated  version  of  the  signal 
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which  when  compared  to  Fig.  2.3.6  shows  how  much  higher  the 
amplitudes  are  and  it  is  impossible  to  see  that  the  signal 
is  periodic.  The  FT  of  the  signal  is  shown  in  Fig.  2.3. 12 
where  it  is  noticed  that  the  spectrum  is  higher  above  the 
zero  db  line  than  in  the  previous  case  because  the  noise 
power  is  large.  Vie  see  that  the  two  peaks  are  separated  by 
a  very  shallow  dip  and  that  there  is  a  peak  on  the  right 
side  so  that  it  looks  like  there  are  three  peaks  in  the 
range  of  frequencies  the  signal  is  known  to  be  in.  The 
solution  to  our  method  gave  coefficients  from  which  the 
spectrum  of  Fig.  2.3*13  is  obtained.  The  normalized  fre¬ 
quencies  of  the  peaks  are  f,  =0.2143  and  f2=0.2304,  Co¬ 
mparing  these  values  with  the  actual  ones  we  see  that  the 
first  peak  occurs  at  the  exact  frequency  and  the  second 
peak  is  only  slightly  in  error.  This  error  is  due  to  the 
TBW  and  noise.  We  see  that  by  increasing  the  noise  five 
fold  we  are  able  to  get  as  good  results  as  when  the  noise 
was  smaller.  The  poles  and  zeros  from  the  coefficients  are 
poles:  0.219+J0.973  zeros:  0. 177+j0.982 
0.128+j0.994  0.213 

0.0 

which  when  compared  to  Fig.  2.3.5  are  slightly  off.  From 
the  location  of  the  poles  and  zeros  we  take  the  inverse  Z 
transform  of  the  parametric  model  to  obtain  the  equation  of 
the  normalized  signal  s(n) 
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s(n)=0.435cos( 1.349n+0. 0609) +0.570cos(1.4478n+0. 0961) 

(2.3-3) 

which  when  compared  to  (2.3-1)  shows  that  the  ratio  of  the 
amplitudes  is  incorrect  and  that  the  solution  introduces 
phases  which  were  not  there  originally.  This  signal  s(n) 
is  shown  in  Fig.  2.3*15  where  the  number  of  periods  of  the 
envelope  of  the  signal  is  eight  compared  to  the  actual 
number  of  six. 

In  the  next  case,  Fig.  2.3*16  shows  the  time  truncated 
version  of  the  signal  for  a  TBW=0.5  and  noise  with  var¬ 
iance  =  1.0.  The  spectrum  in  Fig.  2.3*17  shows  that  the 
first  peak  is  not  distinguishable  because  the  dip  between 
the  two  peaks  is  at  the  same  level  as  the  first  peak. 

With  our  method,  the  spectrum  obtained  is  shown  in  Fig. 

2.3*18.  It  shows  two  peaks  located  at  Ti  =0.2148  and 
f2 =0.2285  which  are  at  the  same  locations  as  in  the  first 
case.  By  decreasing  TBW  our  method  of  solution  was  able  to 
get  as  good  results  as  when  the  TBW  was  larger  and  the  two 
peaks  were  more  clearly  defined.  Again  we  notice  that  the 
first  peak  is  larger  which  is  incorrect  because  of  the 

location  of  the  poles  and  zeros  in  the  z-plane.  The  lo¬ 

cation  of  the  poles  and  zeros  are 

poles:  0. 219+j0.976  zeros:  0. 173+J0.979 
0.135+j0.991  0.056 


which  have  been  plotted  in  Fig  2.3*19*  One  of  the  biggest 
errors  in  this  case  is  the  location  of  the  zero  on  the  real 
axis.  This  is  causing  the  amplitudes  to  be  incorrect.  By 
taking  the  2  transform  of  the  parametric  model  using  these 
poles  and  zeros,  the  time  signal  obtained  is  shown  in  Fig. 
2.3*20.  The  envelope  of  the  signal  has  seven  periods  in 
the  time  shown.  The  normalized  signal  in  Fig.  2.3*20  can 
be  v/ritten  as 

s(n)=0.544cos( 1.3499n-0.C423)+0.469cos( 1 ,4353n+0.0924) 

(2.3-4) 

As  a  final  case  we  chose  to  have  a  TBW=0.359  with  a  noise 
variance  of  cr2=1.0.  The  short  duration  signal  is  shown  in 
Fig.  2.3.21.  The  spectrum  of  the  signal  is  shown  in  Fig. 
2.3*22  where  we  see  very  broad  peaks  throughout  with  one 
very  broad  peak  where  the  two  pulses  should  be.  This  broad 
peak  is  made  up  of  the  two  pulses  smeared  together  because 
of  the  very  small  TBW.  Looking  at  it,  one  could  not  tell 
that  two  pulses  belong  there.  Therefore,  to  see  where 
these  two  peaks  occur  we  used  our  method  to  determine  the 
coefficients  of  the  parametric  model.  The  resultant  spec¬ 
trum  is  shown  in  Fig.  2.3.23.  The  normalized  frequencies 

of  the  pulses  are  f|  =0.2148  and  72=9*2305*  The  second 
estimate  is  incorrect  but  is  the  identical  location  found 


when  the  TBW=0.7l8  and  the  noise  variance  was  a  -5*0.  By 
decreasing  the  TBW,  our  method  did  no  worse  than  for  the 
second  case  investigated.  The  location  of  the  poles  and 
zeros  of  Fig.  2.3*24  are 

poles:  0.219+j0.9?6  zeros:  0. 159+j0*9?8 

0. 122+ jO. 992  0.277 

0.0 

By  taking  the  inverse  Z  transform  of  the  pole- zero  model, 
we  found  the  normalized  estimate  of  the  signal  to  be 


s (n) =0. 608cos ( 1. 3496n+0. 1452 )+0.404cos(1.4483n-0. 1072) 

(2.3-5) 


which  is  plotted  in  Fig.  2.3*25* 


2.4  Conclusion 

We  have  shown  a  method  to  determine  the  location  of  two 
frequencies  when  the  time  bandwidth  product  is  small  and 
when  noise  is  added  to  the  signal  whose  variance  is  equal 
to  or  greater  than  the  power  of  the  smallest  cosine.  The 
method  makes  use  of  the  representation  of  two  impulses  in 
the  z-plane  which  is  used  to  obtain  constraints  for  the 
maximum  and  minimum  values  the  coefficients  of  our  para¬ 
metric  model  can  take.  These  constraints  along  with  the 


relationship  these  coefficients  have  in  the  spectrum,  we 
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used  linear  programming  to  solve  the  set  of  linear  equa¬ 
tions  by  minimizing  the  1,  norm.  Once  the  coefficients 
have  been  solved  for  we  are  able  to  obtain  the  normalized 
estimate  of  the  time  signal  by  inverse  Z  transforming  the 
parametric  model. 

Four  cases  were  used  to  simulate  our  method  in  showing  how 
affective  it  is.  It  was  found  that  by  decreasing  the  TEW 
we  obtained  similar  results  to  when  a  higher  TBW  was  used 
v/ith  a  larger  noise  power.  In  the  last  case  examined  the 
TBW  was  so  low  that  the  two  pulses  in  the  spectrum  smeared 
together  showing  a  very  broad  peak.  Our  method  was  able  to 
obtain  a  good  estimate  of  the  frequencies.  The  time  es¬ 
timate  of  the  signal  is  different  because  the  location  of 
the  poles  and  zeros  is  not  exact. 
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Fig.  2.4.2  Log  magnitude  of  spectrum  of.s(n) 
having  no  noise  and  all  512  points  are  used  in 
evaluating  the  Fourier  transform. 
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Fig.  2.4. 
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Fig.  2.4.7  Log  magnitude  of  spectrum  when  the 
TBW=0.719  and  noise  has  variance  0-2  =  1.0. 
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Fig.  2.4.13  Log  magnitude  of  the  estimated  spe¬ 
ctrum  using  new  method  with  positive  peaks  oc- 
curing  at  the  normalized  frequencies  x| =0.2148 
and  f2 =0.2304. 


*  &o  mo  mo  sqs.0  mo 

nrasRm-ES 


Fig.  2.4.16  Signal^  s(n)  when  T3W=0.5  and  the 
noise  has  variance  cr2  =  i  o. 
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Fig.  2.4.1s  Log  magnitude  of  the  estimated  spe¬ 


ctrum  using  the  new  method  with  positive  peaks 
occuring  at  normalized  frequencies  T|  =0.2148  and 
?*2  =  0. 2285  • 


Fig.  2.4.20  Estimated  signal  sin;  obtamea  T>y 
taking  the  inverse  Z-transform  of  the  transfer 
function  made  up  of  the  estimated  poles  and 
zeros . 


Fig.  2.4. 2i  Short  time  duration  signal  s(n) 
where  TBVI=0.359  and  the  noise  variance  is 

or2  =  1 . 0 . 


Fig.  2.4.22  Log  magnitude  of  spectrum  of 
for  the  short  time  duration  signal  s(n). 


Fig.  2.4.24  Locations  of  the  estimated  poles  and 
zeros  in  the  first  quadrant  of  the  z-plane. 
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Chapter 


FINITE  TIME- BANDWIDTH  PRODUCT 
HILBERT  TRANSFORM 

3.0  Introduction 

When  analyzing  a  process,  one  finds  that  under  certain 
circumstances,  the  real  and  imaginary  parts  of  a  signal  are 
related  through  a  specific  relationship,  or  the  amplitude 
and  phase  are  related  by  that  same  relationship.  In  dif¬ 
ferent  disciplines,  the  relationships  are  known  under 
different  names.  In  mathematics  literature  these  relations 
are  referred  to  as  Poisson's  formulas  [32]  »  in  optics  they 
are  known  as  the  dispersion  relations  [33]  »  and  in  signal 
processing  theory  the  relations  are  called  Hilbert 
transform  relations. 

The  Hilbert  transform  is  used  in  the  phase  retrieval  prob¬ 
lem  which  arises  when  the  wave  phase  is  apparently  lost  or 
impractical  to  measure  and  only  intensity  data  is  avail¬ 
able.  This  situation  occurs,  for  example,  in  electron 
microscopy  where  the  index  of  refraction  structure  of  thin 
films  or  the  height  distribution  of  a  surface  is  to  be 
determined  from  the  intensity  distribution  in  the  far 
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field.  The  phase  problem  also  occurs  in  coherence  theory 
[3^],  signal  processing  [35.36,37].  antenna  array  design 
[38],  Fourier  transform  spectroscopy  [39],  and  design  of 
radar  signals  [4o] . 

If  we  are  given  a  real  signal  f(t)  that  is  square  inte- 
grable,  and  bandlimited  to  a  frequency  Si,  then  the  signal 
f(t)  can  be  represented  as 

f(t)=Re  f*(t)expj  Sit , 

where 


By  definition,  the  Fourier  transform  of  f(t)  is  truncated 
at  w=-Xl.  That  is,  F(w)  is  identically  zero  for  w<-&. 
Multiplying  f(t)  by  exp  j  SI  t  shifts  the  spectrum  to  the 
right  by  an  amount  SI.  The  resulting  complex  signal  has  no 
negative  frequency  components. 


The  preenvelope  of  a  real  signal  f(t),  is  the  complex- valued 
function 


m(t)=f(t)  +  j?'(t)  (3.0-1) 

The  real  signal  is  the  real  part  of  the  preenvelope  m(t). 
The  preenvelope  is  also  called  the  analytic  signal.  An 
analytic  signal  has  the  property  that  the  envelope  of  f(t) 
is  the  absolute  value  |m(t)|  of  its  preenvelope  which  is  of 
use  in  modulation  theory  [35]* 

In  (3.0-1)  we  have  that  the  real  and  imaginary  parts  of 
the  signal  m(t)  are  related  by  the  Hilbert  transform 
given  by 

f(t)  =  l  P  /•C°f(r)  dr  =  H [f ( t ) ]  (3-0-2) 

tr  -loo  t-  t 


where  P  denotes  the  Cauchy  principle  value  given  by 


f( 


t)  =  l  lim  f  r~f(r)  dT  +  fhr)  dr  1 
T T  lv«t-T  -cot-r  J 


(3.0-3) 


The  frequency  spectrum  of  f(t)  is  given  by  its  Fourier 


transform 


F(w)  =  F[_f(t)J  =  1  f  (t)exp(- j27rft)dt 

•-CD 

(3.0-4) 

and  F(w)  exhibits  real-even,  imaginary-odd 

cause  f(t)  is  real.  The  spectra  of  f(t)  and 

lated  through  the  relationship 

symmetry  be- 

f(t)  are  re- 

F[f(t)]=  -  jsgn(w)F(w) 

(3-0-5) 

where 

’  +  1  ,  w£0 

sgn(w)=-  0  ,  w=0 

-  1  ,  ws 0 

(3.0-6)' 

From  this  the  Hilbert  transformation  of  f(t)  can  be  viewed 

as  f(t)  passed  through  a  -90°  phase-shift  network  whose 

frequency  and  impulse  responses  are 

G(w)=  -jsgn(w) 

g(t)=  l/7Tt 

(3.06-7a) 

( 3. 0-?b) 

Because  G(w)=  -1,  we  must  have 

f(t)=  -HjT(t)] 

(3.0-8) 

and  f(t),  'f(t)  are  termed  a  Hilbert  pair. 
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The  spectrum  of  m(t)  is  shown  to  have  zero  negative  fre¬ 
quency  components  by  substituting  (3*0-5)  into  the  Fourier 
transform  of  (3*0-1).  This  gives  us 

M(w)=F(w)  +  jF[f(t)]  (3*0-9) 

or 

M(w)=  F(w)+j (-jsgn(w)F(w) ) 

=F(w)  [l+sgn(w)]  (3.0-10) 

For  w  less  than  0,  from  (3*0-6)  we  see  that  the  sign  fun¬ 
ction  is  -1  and  the  function  M(w)  disappears.  For  w 
greater  than  0,  M(w)=2F(w),  and  for  w=0  M(w)=F(0).  This 

can  be  rewritten  as 


2F(w)  ,  w>0 
M ( W )  =  *  F(0)  ,  w=0 

0  ,  w<0 


(3*0-11) 


When  the  spectrum  of  M(w)  is  banalimited  to  +Ji,M(w)  van¬ 
ishes  for  w>Xi,  and  the  finite  energy  Em  of  the  signal  is, 
by  Parseval's  theorem 


(3*0-12) 


Let  z=T+j  <j  define  a  complete  variable.  By  inverse 

Fourier  transformation 


m(z)  =  1  J  I 

27 r 

■-/: 


M(w)expjwt  dw 


I>i(w )  exp-w  cr  expjv/T  dw. 


(3* 0-1 3a) 


( 3* 0-1 3b) 


It  becomes  evident  that  given  the  integrability  and  con¬ 
vergence  properties  for  the  existence  of  the  Fourier  trans¬ 
form  that  (3-0-13)  must  converge  for  any  a>0.  Therefore, 
m(z)  must  be  free  of  zeros  in  the  upper  half  of  the 
z-plane . 


If  we  let  m(z)  be  bandlimited  to  £1 ,  from  (3. 0-1 3)  we  have, 


hn(z)i*  *  1  f  M(w)  exp-w^-  dw  1 


(3-0-14) 


From  the  Schwartz  inequality  (3-0-14)  becomes 


|  m  ( z )  |*  -  1  f  M(w)  dw  1  fexp-2wcr  dw 


■-£ 


-  Err\  l-eXp(-2  £1  c T  ) 


4  7T0- 


(3-0-15) 


•. 
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where  Em  is  the  finite  energy  given  by  (3*0-12). 

This  shows  that  m(z)  is  bounded  everywhere  in  the  finite 
lower  half  and  whole  upper  half  of  the  z-plane.  When  M(z) 
is  bandlimited ,  then  m(z)  must  be  an  entire  function.  This 
relationship  does  not  depend  on  bandv/idth  limitations. 
With  n,  =ao,  (3.O-I5)  is  finite  for  any  cr  >0  provided  that 
m(t)  is  square  integrable.  A  finite  bandwidth  insures  that 
the  lower  half  plane  singularities  are  at  infinity.  The 
Hilbert  transform  can  also  be  used  to  describe  the  re¬ 
lationships  between  the  modulus  |n(t)j  and  the  phase  <f>  (t) 
of  an  analytic  signal.  It  is  assumed  that  M(w)  is  band- 
limited.  Equation  (3.0-2)  can  be  written  in  phasor  nota¬ 
tion  as 

m(t)=  m(t)  exp<£(t)  (3*0-16) 

where 

|  m ( t )j  =*v/f2(t)+f2(t)  (4 . 0-  17a) 

(t)=  arctan  f (t)/f (t)  (3.0-17b) 

By  taking  the  complex  logarithm  of  ( 3 . 0— 1 6 )  we  have 


In  m(t)  =  In  | m ( t ) J  +  j<£(t) 


(3*0-18) 


in  which  the  logarithmic  modulus  and  the  phase  are  the  real 
and  imaginary  parts  of  a  time  function.  Certain  conditions 
have  to  be  imposed  in  order  for  the  Hilbert  relation  to 
hold  between  the  modulus  and  phase.  We  first  note  that 
In  m(t)  is  not  square  integrable,  and  in  general,  m(t)  —  0 
as  t  —cd  and  thus  In  m(t)  —  co  as  t  —  co  .  However,  in 
certain  instances,  we  can  modify  m(t)  by  the  addition  of  a 
constant  unit  amplitude  so  that  the  modified  function 
ln(l+m(t))  is  square  integrable  if  |m(t)j  is  chosen  to  be 
less  than  one  [4l] . 

Another  way  to  avoid  the  problem  is  to  study 

ln'm(z)=  d  In  m(z)  (3.0-19a) 

dz 

=  ln*  m(z)  +  j  <£  •  (t)  (3.0-19b) 

=m'(z)/m(z)  (3*0-19c) 

From  the  Payley-Wiener  theorem  [42],  if  m(z)  is  entire  and 
square  integrable  along  the  real  axis,  and  if  M(w)  is  band- 
limited  to  £lrthen 

m(z)=0(exp(  & | z| ))  (3.0-20) 

where  0  means  "order  of"  defined  as  the  maximum  absolute 
value.  This  means  that  the  maximum  absolute  value  of  m(z) 


& 
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is  increasing  exponentially  with  radius  r=jz|.  From  this 
we  have  that  In  m(z)  at  worst  will  vary  linearly  with  |z|, 


and  ln'm(t)  must  therefore  vanish  or  tend  to  a  finite  con¬ 


stant  in  the  extremities  of  the  z-plane. 


Because  m(z)  is  analytic,  ln'm(z)  will  have  upper  half  plane 


singularities  only  if  m(z)  has  zeros  in  the  upper  half 


plane.  If  ln'm(z)  is  square  integrable,  we  then  have  the 


relationships 


4>  '  (t)=  K  [  In'  m(t  )J 
In'  Jm(t)|  =  -H  [  <£  ' (t)] 


( 3. 0-21a) 


( 3 • 0-2  lb ) 


where  m(z)  must  be  zero  free  in  the  upper  half  of  the 


z-plane.  The  relationship  where  the  phase  and  magnitude 


are  related  through  the  Hilbert  transform  is  known  as  the 


minimum  phase  condition  where  the  phase  displacement  is  the 


smallest  possible  for  its  gain. 


From  (3.0-21)  we  can  always  calculate  (t)  and  In  m(t)  by 


integration.  The  results  obtained  will  generally  not  be 


unique.  In  certain  applications  the  fact  that  the  solution 


is  not  unique  does  not  play  a  role  in  the  solution  [35] 


The  theory  just  described  has  been  for  analytic  time  sig¬ 
nals  whose  spectrum  has  zero  negative  frequency 


fey  liC* Lt'. v/'oj.  * 


components.  A  dual  of  this  is  a  time  signal  that  has  no 
negative  time  components.  It  can  be  shown  that  the  real 
and  imaginary  parts  of  the  spectrum  are  related  by  the 
Hilbert  transform.  Signals  that  have  no  negative  time 
components  are  termed  causal.  In  cybernetics,  a  causal 
signal  that  is  also  square  integrable  is  called  a  v/avelet 

[37]  • 

The  solution  of  the  phase  retrieval  problem  can  either  be 
solved  for  by  relying  on  the  analyticity  of  the  signal 
where  the  Hilbert  transform  can  be  used,  or  a  computational 
procedure.  In  the  case  of  what  in  optics  is  called  a  weak 
object  [43,44]  one  has  a  signal  F(w)  given  by 

F(w)  =  1+M(w)  (3.0-22) 

v/here 

J  M(w)|«  1. 


With  this  condition  we  can  approximate  the  real  part  of  the 
signal  M(w)  as 

Re  [M(w)]  «0.5(F(w)-l)  (3-0-23) 

and  the  imaginary  part  of  the  signal  M(w)  as 


■V*. 
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Im  [M(w)J«  hr  [Re  M(w)] 


(3-0-24) 


As  we  can  see  we  have  to  have  the  special  condition  that 


the  signal  M(w)  has  a  constant  dc  offset  and  the  signal 


magnitude  is  much  smaller  than  the  offset. 


Another  procedure  is  the  apodizing  technique.  In  this 


technique,  the  function  m(t)  is  modified  so  that  the  zeros 


of  the  modified  function  are  displaced  so  that  their  con¬ 


tribution  to  the  phase  is  diminished.  The  logarithmic 


Hilbert  transform  can  then  be  used  to  calculate  the  chase. 


The  apodization  concept  follows  from  considerations  of  the 
bandlimited  Fourier  transform  [33]: 


m  & 

z)=  J  f(t)  exp(-yt)cos(xt+arg  f(t))dt 


■  j  Jo  f(t)  exp(-yt )sin(xt+arg  f(t))dt 


(3.0-25) 


where  z=x+jy.  If  f(t)  is  made  to  decrease  more  rapidly, 
near  the  lower  limit  of  the  integrals  then  larger  x,y  must 
be  required.  This  means  that  a  larger  zero-free  area  about 
the  origin  is  created,  by  reducing  the  phase  contributions 
of  the  zeros  in  the  lattice  by  having  f(t)  decrease  more 


rapidly  near  the  edges  of  the  window.  This  c.-.n  be  accom¬ 


plished  by  multiplying  f(t)  by  a  suitable  filter. 


Nakajima  and  Asakura  [45]  have  used  the  Gaussian,  sine  and 
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triangular  filters  to  perforin  the  apodization  which  is  a 
multiplication  of  the  autocorrelation  of  the  signal  with 
the  desired  filter  function.  The  autocorrelation  is  the 
inverse  Fourier  transform  of  the  magnitude  in  the  fre¬ 
quency  domain. 

As  was  already  mentioned,  it  becomes  difficult  to  obtain 
the  phase  from  magnitude  measurements  if  the  function  has 
zeros  in  the  upper  half  of  the  s-plane.  Procedures  have 
been  developed  to  determine  where  in  the  upper  half  plane 
the  zeros  occur  so  that  they  could  be  flipped  into  the 
lower  half  plane.  By  accomplishing  this  flip,  the  mag¬ 
nitude  in  the  frequency  domain  remains  the  same.  Nakajima 
and  Asakura  [46]  have  devised  a  method  to  determine  the 
position  of  the  zeros  from  the  magnitude  of  the  signal  and 
the  magnitude  of  the  Fourier  transform  of  the  signal. 
From  this  along  with  the  logarithmic  Hilbert  transform  and 
a  nonlinear  least-squares  parameter  estimation  technique, 
they  have  been  able  to  determine  the  phase. 

Gerchberg  and  Saxton  [47,48]  were  first  to  suggest  the 
use  of  both  the  magnitude  of  the  signal  and  the  magnitude 
of  the  spectrum  to  obtain  the  phase.  Misell  [  49]  proposed 
the  use  of  the  intensity  distributions  of  two  slightly 
defocused  images.  In  the  Gerchberg-Saxton  algorithm,  at 
each  step  the  computed  values  of  the  intensity  are  cor- 


rected  by  combining  with  the  measured  values.  Other  sug¬ 
gestions  [50]  for  use  of  the  magnitudes  of  the  signal  and 
the  spectrum  have  involved  the  solution  of  algebraic  equa¬ 
tions  which  connect  the  sampled  magnitudes  in  time  and 
frequency  with  coefficients  of  the  discrete  Fourier  trans¬ 
forms  of  the  function. 

Bates  et  al.,  [ 51 • 52, 53]  have  obtained  a  simple  algebraic 
derivation  to  the  phase  problem,  and  have  devised  an  easily 
implementable  algorithm.  The  algorithm  works  for  images 
that  are  weakly  localized.  What  this  means  is  that  the 
energy  of  each  sample  is  slightly  spread  about  the  sample. 

An  alternative  approach  to  phase  retrieval  problem  is  to 
employ  one  of  the  gradient  search  methods.  It  has  been 
shown,  [5^] ,  that  one  such  method,  the  steepest-descent 
method,  is  closely  related  to  the  error  reduction  algo¬ 
rithm.  Its  one  drawback  is  that  it  converges  slower  than 
other  gradient  search  methods  that  are  available. 

Much  work  has  recently  been  done  in  reconstructing  the 
signal  from  its  phase  [60,61,62,63].  This  is  because  it  is 
possible  to  relate  the  phase  of  a  signal  to  the  signal 
under  certain  conditions,  which  makes  it  easier  to  perform 
than  to  obtain  the  signal  from  its  magnitude. 


In  this  chapter  we  describe  a  new  form  of  the  Hilbert  trans¬ 
form  which  is  based  on  finite  time  observations  of  a 
signal.  It  will  be  shown  that  the  kernel  does  not  have  the 
inherent  singularity  that  the  regular  Hilbert  transform 
exhibits.  l.'e  will  show  that  the  modified  Hilbert  trans¬ 
form  has  similar  properties  of  the  regular  Hilbert  trans¬ 
form.  The  modified  Hilbert  transform  will  be  used  with  a 
method  that  we  call  the  method  of  partitioning  to  de¬ 
termine  the  phase  of  a  function  given  its  magnitude.  Nu¬ 
merical  results  will  be  used  to  show  how  this  new  method 
works . 
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where  R(w)  is  the  real  part  and  X(w)  is  the  imaginary  part 
of  the  complex  signal.  By  looking  at  ( 3 • 1—1 )  and  (3.1-2) 
we  see  that  the  integrals  are  evaluated  in  both  the  time 
and  frequency  domains,  and  because  of  the  even  symmetry  the 
integrals  are  also  evaluated  for  only  positive  values. 

If  we  were  to  assume  that  the  signal  is  observed  in  the 
time  domain  for  a  finite  time  T,  then  in  (3*1-1)  and 
(3.1-2),  the  upper  limit  of  integration  for  t  is  T.  For 
simplicity  we  will  obtain  the  modified  Hilbert  transform 
from  (3*1-2)  with  the  understanding  that  the  procedure  is 
applicable  to  (3. 1-1).  With  T  as  the  upper  limit,  (3.1-2) 
can  be  written  as 


;v 


CO 


X(w)=-2  f  R(y)  f  sinwt  cosyt  dt  dy 

*  J 0  Jo 


(3*1-3) 


We  see  from  (3*1-3)  that  the  integral  evaluated  over  t  can 
be  evaluated  to  give  the  equation  only  in  terms  of  w  and  y. 
The  integration  gives  us  the  following 


l  f  si 

Jo 


2  f  sinwt  cosyt  dt=  l-cos(w+y)T  -  l-cos(w-y)T  (3*1-4) 


(w+y) 


(w-y) 


which  can  be  placed  back  into  (3*1-3)*  Therefore,  we  have 
from  (3*1-3)  and  (3.1-4) 
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x(w)=-l  r  R(y ) 


CO 


l-cos(w+y)T  +  l-cos(w-y)T 


(w+y) 


(w-y) 


dy 


(3.1-5) 


If  the  complex  signal  is  assumed  to  be  bandlinited  to  w=  Si 
we  can  rewrite  (3.1-5)  as 

X(w)=-1_  r  R(y)  l-cos(w+y)T  dy 

*  Jo 


(w+y) 


-1_  f  R(y )  l-cos(w-y)T  dy 

^ T  Jo  - 

(w-y) 


(3.1-6) 


where  we  have  broken  the  integral  into  the  sum  of  two  inte¬ 
grals.  The  two  integrals  can  be  combined  into  one  by 
changing  the  variable  in  the  first  integral  to  give  the 
final  result 


X(w)=-JL  f  R(y)  l-cos(w-y)T  dy 

W  J-A 


(3.1-7) 


(w-y) 


The  above  equation  can  also  be  written  with  different 
limits  of  integration.  It  is  seen  that  both  the  variables 
w  and  y  take  on  the  values 

-SI  =w=fl 
-  fl=y=  SI. 
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If  new  variables  are  defined  given  by 

w '  =w/fl 
y  ,=y/il 

we  can  rewrite  the  inequalities  as 

— l=w’=l 

-l=y'=l 

With  these  new  variables  (3* 1-7)  can  be  written  in  the  form 


X(flW)-M  f  R(fly')  1-cos ftT(w'-y' )  dy' 

(w  *  -y ' ) 


(3.1-8) 


A  reason  why  one  may  be  interested  in  this  form  of  the 
equation  is  that  the  time-bandwidth  product  is  directly 
incorporated  into  the  equation.  This  is  the  form  of  the 
transform  we  have  termed  the  modified  Hilbert  transform. 
From  (3.1-1)  we  can  go  through  the  same  analysis  to  obtain 


R(w)  =-_l_  |  X(y)  l-cos(w-y)T  dy 

(w-y) 

R(nw)-if  X(^ly')  1-cos  ft  T(w'-y')  dy' 

wL 


(3.1-9) 


(3. l-9b) 


(w'-y’ ) 
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Equations  (3*1-7)  and  (3*l-9a)  can  be  rewritten  using  tri- 
ginometric  identities  as 


x(w)=-l_  j 

w)=-l_  f* 
irj-i 


R(w)  = 


'R(y)  sin20.5T(w-y)  dy 
i  0 . 5 (w-y ) 

X(y)  sin20.5T(w-y)  dy 
1  0.5 (w-y) 


(3. l-10a) 


( 3. 1-1  Ob) 


which  if  we  take  a  closer  look  at  the  kernel,  contains  the 
sampling  function  given  by 


sinc0.5Tw  =  sinO . 5Tw/( 0 . 5Tw) 


(3*1-11) 


This  sampling  function  is  a  consequence  of  the  finite  time 
observation  of  the  signal.  We  see  that  the  modified 
Hilbert  transform  is  a  convolution  of  a  signal  with  a  ker¬ 
nel  evaluated  over  a  finite  bandwidth.  The  difference 
between  the  regular  and  modified  Hilbert  transform  is  the 
kernel  given  by 


K(w)  =  1-cosTw  =  TsinO.5wTsincO.5wT 


(3*1-12) 


which  differs  with  the  kernel  of  the  Hilbert  transform 
given  by 
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K ’ (w) =  1/w.  (3.1-13) 

In  (3*1-13)  we  see  that  the  kernel  has  a  singularity  at  w=0 
where  the  function  approaches  infinity.  For  the  modified 
Hilbert  transform  kernel  we  see  that  in  the  limit  as  w 
approaches  to  zero,  the  sampling  function  approaches  one 
and  the  sine  approaches  zero  giving  the  overall  limit  as 
zero.  The  kernel  of  the  modified  Hilbert  transform  has 
an  envelope  which  is  identical  to  the  kernel  K'(w)  with 
sinusoidal  variations  which  makes  it  go  to  zero  as  w  goes 
to  zero.  In  Fig.  3*1.1  we  see  the  kernel  K(w)  for  T=1  and 
f=12  where  we  have  that  w=2irf.  We  see  that  for  positive 
frequencies  the  kernel  has  only  positive  values  and  for 
negative  frequencies  only  negative  values  as  in  the  case 
of  K'(w).  In  Fig.  3*1.2  we  have  expanded  a  portion  of  the 
spectrum  and  evaluated  the  kernel  for  frequencies  up  to  f=3 
with  T=l.  From  this  and  the  previous  figure  we  note  that 
the  majority  of  the  energy  is  located  in  the  lower  fre¬ 
quencies. 

By  a  change  in  variables  we  can  get  the  equivalent  forms  of 
(3.1-10a)  as 

X(w)=-_1_  /*nsinc0.5(w+y)T  sin0.5(w+y)T  R(y)  dy  (3.1-l4a) 

"J-SI 

X(w)=-^  I  sincO.5ytsinO.5yT  R(w-y)  dy 

*  J-n 


v.v.v.a  ,vy/v.  v.v.- .v.-  .y. 


(3. l-l4b) 
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which  can  be  extended  to  (3.1- 10b). 


To  evaluate  the  integral  we  replace  the  integral  by  a  sum¬ 
mation  to  give  us  (3.1- 10a)  in  the  form 


M+l 


X(A,  )=-l  V  (l-cos((Aj  -  A;  )T)  R(A.  )A 

Aj  -  As 


(3.1-15) 


7T 


i*l 


where  M  is  the  number  of  intervals  the  integral  is  broken 
into,  and 


A  =  zCl/m 
Aj  =  ( j-l)A  -CL 

Ai  =(i-l )A  -  Cl  i,  j  =  l,2 . M 


(3.1-16) 


Using  (4.1-16)  we  can  rewrite  (3.1-14)  as 


M+l 


X(Aj  )=-J_  y  ( l-cos(  ( j-i)A  T) )  R(A;) 

( j-i) 


(3.1-17) 


To  rewrite  (3. 1-1 Ob)  the  X  and  R  need  only  be  interchanged. 


As  an  example  to  show  how  this  modified  Hilbert  transform 
evaluates  the  imaginary  part  from  the  real  part  of  the 
Fourier  transform  of  a  causal  signal  we  have  chosen  our 
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signal  f(t)  to  be  a  unit  amplitude  step  for  a  duration  T. 
Therefore,  we  have  f(t)  given  as 

f(t)=  f  1  , 0=t=T 

-  lo  , otherwise  (3. 1-18) 

The  Fourier  transform  of  f(t)  is  given  by 

F(w)=0.5Tsinc0.5wTexp(- jO . 5wT) 

=0.5TsincC.5wT(cos0.5wT- jsin0.5wT)  (3* 1—19) 

Therefore,  if  R(w)  is  the  real  part  c “  the  signal  given  by 

R(w)=  0 . 5TcosO . 5wTsincO . 5wT  (3.1-20) 

then  the  modified  Hilbert  transform  should  yield  the  im¬ 
aginary  part  X(w)  given  by 

X(w)=  -0.5Tsin0.5wTsinc0.5wT  (3.1.21) 

In  Fig.  3.1-3*  the  real  part  of  the  signal  R(w)  is  given 
by  the  solid  line.  It  is  seen  to  be  an  even  function.  To 
obtain  the  imaginary  part  X(w),  (3. 1-17)  was  used  with 

M=256,  T=  1  and  f=£l/27r=12.  The  imaginary  part  is  drawn 
using  a  dashed  line.  To  see  how  closely  the  result  cor¬ 
responds  to  the  actual  values,  (3. 1-21)  was  also  drawn 
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using  a  dotted  line.  We  see  that  the  actual  and  the  cal¬ 
culated  imaginary  part  are  overlapping,  indicating  that  the 
calculated  function  is  identical  to  the  actual  one.  As  a 
convenience  for  graphing,  all  of  the  functions  will  be 
normalizing.  The  normalized  constant  is  obtained  by  find¬ 
ing  the  maximum  «.  olute  value  of  the  function  and  dividing 
the  function  by  is  value.  This  assures  us  that  the  mag¬ 
nitude  of  the  function  will  always  be  less  than  or  equal  to 
one.  In  Fig.  J.l.J,  as  will  always  be  the  case,  the  nor¬ 
malizing  constant  KR  for  the  real  part  of  the  signal  is 
given  by  0.5*  The  normalizing  constant  for  the  actual 
imaginary  part  KI  is  given  by  0.J62,  and  the  normalizing 
constant  of  the  calculated  or  estimated  function,  KE  is 
given  by  O.360.  From  the  normalizing  constants  KI  and  KE 
we  see  again  that  the  results  obtained  using  the  modified 
Hilbert  transform  are  very  close. 

In  (3.I-I6)  we  have  that  the  spacing  between  samples  is 
A  =  2&/M  =  4ir  f/M  (3.  l-22a) 


or 


A'  =  A/2  ir  =  2  f/M 


(3. l-22b) 


From  the  sampling  theorem  for  the  Hilbert  transform  T5?]  » 
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we  have  that 


A '  =  l/T 


(3.1-23)) 


from  which  we  can  obtain  an  expression  for  the 
tine-bandwidth  product  TBP  associated  with  the  sampling 
interval  A'.  From  (3.1-22)  we  can  deduce  that  the  band¬ 
width  B  is  given  by 


B  =  mA’/2  =f 


From  (3*1-23)  and  (3*1-24)  we  have 


f  =  M/2T 


(3.1-24) 


(3*l-25a) 


fT  =  M/2 


(3.l-25b) 


Equation  (3.1-25b)  states  that  if  one  is  to  sample  accor¬ 
ding  to  (3.I-23)  then  the  time-bandwidth  product  is  equal 
to  M/2.  In  the  case  where  the  time-bandwidth  product  is 
less  than  M/2,  the  signal  is  being  oversampled,  and  if  the 
product  is  greater  than  M/2,  the  signal  is  undersampled. 
In  speech  analysis  [58,59]  »  for  certain  processing 
techniques,  one  finds  it  necessary  to  oversample  the  signal. 
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In  the  example  of  Fig.  3.1.3,  the  TBP=12  which  tells  us 
that  we  are  oversampling  by  10. 67  times.  To  test  the  modi¬ 
fied  Hilbert  transform  in  how  well  it  works  for  calcu¬ 
lating  the  real  part  of  the  signal  from  the  imaginary  part 
with  identical  constants  as  before,  we  use  the  negative  of 
(3.1-21)  to  see  how  closely  our  results  match  (3.1-20). 
Fig.  3.1.4  shows  the  results  where  we  see  that  the  graph  of 
the  actual  real  part  of  the  signal  given  by  the  dotted  line 
and  the  graph  of  the  calculated  real  part  are  overlapping. 
The  imaginary  part  of  the  signal  is  given  by  the  solid 
line,  where  for  quick  identification  we  note  that  the  imag¬ 
inary  part  is  odd.  The  normalizing  constants  for  Fig. 
3.1.4  are  KI=0.326,  KR=0.5  and  KE  the  normalizing  constant 
of  the  estimate  to  the  real  part  of  the  signal  is  0.494. 
With  the  shape  of  the  estimate  being  identical  to  the  shape 
of  the  actual  signal,  and  having  the  normalizing  constants 
in  close  agreement,  we  can  say  that  the  estimate  is  exact. 


In  the  above  examples  the  bandwidth  of  the  signal  was  taken 
to  be  12.  This  bandwidth  gave  us  the  limits  of  inte¬ 
gration,  and  we  observed  that  the  calculated  estimates  of 
the  real  and  imaginary  parts  of  the  signal  were  almost 
identical  to  the  actual  values.  To  see  what  effect  re¬ 
ducing  the  bandwidth  would  have  on  our  solution,  we  next 
chose  to  have  the  bandwidth  f=3*  Again  we  have  M=256  and 
T= 1 .  The  TBP  here  is  given  by  3  so  that  we  are  oversampling 
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by  42.67  times.  Equation  (3.1-20)  was  used  as  the  input 
and  we  wanted  to  see  how  closely  our  results  using  the 
modified  Hilbert  transform  would  be  to  (3*2-21) .  In  Fig. 
3.1.5  the  solid  line  is  the  function  given  in  (3.1-20) 
with  KR=0.5*  The  dotted  line  is  the  actual  solution  which 
is  given  by  (3*1-21),  and  its  normalizing  constant  is 
KI=0.362.  The  solution  obtained  by  the  modified  Hilbert 
transform  is  shown  by  the  dashed  line  with  the  normalizing 
constant  given  by  0.356.  We  notice  that  at  the  edges  of 
our  window,  the  estimate  does  not  correspond  to  the  actual 
values.  This  phenomenon  is  caused  by  the  windowing.  The 
real  part  of  the  signal  does  not  have  all  of  the  energy  or 
most  of  the  energy  concentrated  in  the  region  up  to  f=3* 
The  same  can  be  said  for  the  kernel  of  the  transform  which 
we  have  shown  in  Fig.  3*1*2.  By  taking  the  convolution  of 
this  signal  with  the  kernel  over  such  a  small  region,  the 
effect  is  to  produce  smearing  at  the  limits  of  inte¬ 
gration.  The  esimate  about  the  origin  is  seen  to  be  very 
close  to  the  actual  values  because  in  this  region  the  con¬ 
volution  involves  the  main  lobes  of  the  functions  which 
have  enough  energy  to  produce  good  estimates.  In  this 
situation,  we  notice  that  the  sampling  interval  or  con¬ 
sequently  the  TBP  did  not  have  an  effect  on  our  results. 
It  was  our  choice  of  bandwidth  that  caused  the  smearing  at 
the  edges. 
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The  next  example  that  we  wanted  to  choose  was  a  causal 
signal  that-  would  have  all  of  its  energy  in  the  short  time 
that  it  would  be  observed.  For  this  we  chose 


f(t)  =  12exp(-6t )uT(t ) 


(3.1-26) 


where 


u-r(t)  =  f  1  ,  O-t-T 


’  1  . 

« 

.  0  , 


elsewhere 


The  Fourier  transform  of  (3.1-26)  is  calculated  to  be 


F ( w )  =  12(1  -  exp-(6+jw)T) 

6+ jw 


(3-1-27) 


By  choosing  T=l,  the  magnitude  of  the  exponential  is 
2.48xl0-3  so  that  (3*1-27)  can  be  approximated  by 


F(w)  = 


6+ jw 


36+ws 


-j  12w 
36+wJ 


(3.1-28) 


With  this  choice  of  T,  we  find  that  99*9^  of  the  energy  is 
in  the  signal.  This  satisfies  the  condition  that  we  wanted 
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for  our  signal,  in  that  most  of  the  energy  is  contained 
in  the  signal  for  the  short  time  that  it  is  observed.  Next 
we  wanted  to  calculate  the  bandwidth  so  that  95%>  of  the 
energy  be  contained  up  to  that  frequency.  The  energy  in 
the  frequency  domain  is  given  by 

E  =  (2/ 7r  )arctan(27r  f/6)  (3.1-29) 

The  signal  in  (3*1-26)  has  been  normalized  so  that  its 
total  energy  equals  one.  Therefore,  from  (3.1-29)  we  find 
that  for  f=12,  9k. 9%  of  the  energy  is  contained  in  that 
part  of  the  spectrum.  Using  the  real  part  of  (3-1-28)  we 
wanted  to  see  how  well  the  calculated  imaginary  part  is  to 
the  actual  values.  With  M=256,  f=12  and  T=l,  Fig.  3. 1.6 
shows  the  results  obtained.  The  solid  line  is  the  real 
part  of  (3*1-28)  with  a  normalizing  constant  KR=2.0.  The 
actual  imaginary  part  shown  by  a  dashed  line,  overlaps  to 
look  as  if  the  graph  is  one  solid  line.  The  normalizing 
constants  are  KI=1.0  for  the  actual  plot  and  KE=1.0  for  the 
estimate  plot.  We  see  that  the  modified  Hilbert  transform 
gives  exact  results  if  the  limits  of  integration  are  chosen 
such  that  most  of  the  energy  of  the  signal  is  located  be¬ 
tween  the  limits.  From  this  example  we  can  assume  that 
the  bandwidth  can  be  chosen  so  that  95 %  of  the  energy  is 
located  in  that  region.  This  way  the  integration  does  not 
have  to  be  done  over  the  entire  frequency  spectrum. 


In  the  next  section  we  will  use  the  logarithmic  form  of  the 
modified  Hilbert  transform  to  evaluate  the  phase  of  a 
signal  from  its  magnitude.  We  will  also  show  a  new  method 
of  partitioning  the  spectrum  to  obtain  better  estimates  of 
the  phase. 

3.2  Phase  Retrieval  Problem 

In  section  we  discussed  how  for  a  signal  m(t)  there 

exists  a  relationship  between  the  real  and  imaginary  parts 
of  the  complex  logarithm  of  m(t).  If  m(t)  is  bandlimited, 
and  it  can  be  written  in  phasor  notation  as 

m(t)=  Jm(t)J  expj  (t)  (3.2-1) 

then  by  taking  the  complex  logarithm  of  (3.2-1)  we  have 

In  m(t)=  In j m ( t ) j  +j<£(t)  (3.2-2) 

The  ln|m(t)|  is  the  real  part  of  the  signal  and  <j>  (t)  is  the 
imaginary  part.  From  this  we  then  have  that  the  phase 
is  related  to  the  log  of  the  magnitude  by  the  Hilbert  trans¬ 
form.  In  the  notation  of  the  modified  Hilbert  transform  we 


have 


.98 


4>  (w)  = 


-_1_  rftln)A(y)J 
•'-n 


l-cos(w-y)T  dy 


w-y 


(3.2-3) 


The  above  equation  relates  the  phase  of  the  spectrum  with 
the  magnitude  of  the  spectrum  when  the  time  signal  is 
causal  and  real.  If  we  are  dealing  with  an  analytic  sig¬ 
nal,  where  the  spectrum  is  real  and  contains  only  positive 
frequency  components,  then  the  logarithmic  modified 
Hilbert  transform  is  given  by 


t)=  -_1_  JTln|m(r)|  1 


-cos(t-r)  dr 


(3.2-4) 


t-  T 


As  an  example  we  will  use  the  one-sided  exponential  pulse 
given  by  (3.1-22).  As  a  matter  of  convenience,  we  rewrite 
the  Fourier  transform  of  this  signal  when  T=1  as 

F(w)=  12  (3.2-5) 

6+ jw 

The  magnitude  and  the  phase  of  (3.2-5)  can  be  calculated 
to  give 


I.-  * 


|a(w)|  =12/(36+w8  )°9 
^(w)=  arctan(-w/6) 


(3.2-6a) 

(3.2-6b) 
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From  (3.2-6a)  we  notice  that  the  magnitude  approaches  zero 
when  w  approaches  infinity.  This  tells  us  that  the  log  of 
the  amplitude  is  not  square  integrable  as  is  the  necessary 
condition  for  using  the  logarithmic  form  of  the  Hilbert 
transform  in  determining  the  phase.  In  the  limit  as  w  goes 
to  infinity  we  have  that  the  log  of  the  magnitude  also 
goes  to  infinity,  and  we  get 


f  } In  1 A(w) J  jz  dw  — » 
J-  CD 


(3*2-7) 


In  the  previous  section  we  found  that  if  f  is  chosen  to  be 
12,  that  9^.9^  of  the  energy  is  included  in  the  spectrum. 
It  is  with  these  limits  that  we  would  like  to  determine  the 
phase  when  using  the  logarithmic  modified  Hilbert  trans¬ 
form.  If  we  take  the  bandwidth  of  the  signal  to  be  12  we 
have  that  with  =2  it f 


[*)  1«JA(w] 

J-n. 


dw  <  c» 


(3-2-8) 


because  within  these  limits  the  log  of  the  magnitude  never 
goes  to  infinity.  By  this  we  have  bandlimited  the  signal. 
In  Fig.  3*2.1,  the  solid  line  represents  the  log  mag¬ 
nitude,  with  a  normalizing  constant  given  by  KR=1.84l.  The 
dotted  line  is  the  phase  given  by  (3.2-6b)  with  a  nor- 
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malizing  constant  KIM. 491.  The  estimate  of  the  phase 
using  the  logarithmic  modified  Hilbert  transform  is  shown 
with  the  dashed  line.  Its  normalizing  constant  is  given  by 
KE=2.56l.  From  this  figure  we  see  that  the  estimated 
phase  is  incorrect.  There  are  two  reasons  why  the  estimate 
is  incorrect.  The  first,  which  is  hardly  noticeable  is 
that  we  have  bandlimited  the  function  so  that  there  may  be 
some  smearing  of  the  estimate.  The  second  reason  which  is 
the  major  one  is  that  the  magnitude  has  a  pole  in  the  upper 
half  of  the  z-plane.  This  can  be  seen  if  the  log  of  the 
magnitude  is  written  as 

In  A(w)  =0.5  In  (l44/(36+w2  )) 

=0.5  In  ( 144/(6+ jw) (6- jw) )  (3.2-9) 

From  this  we  see  that  a  pole  occurs  at  w=j6  which  is  in  the 
upper  half  of  the  z-plane.  It  is  because  of  this  pole 
that  the  function  is  not  square  integrable  and  leads  to 
such  poor  results. 

In  Fig.  3*2.2  we  have  decided  to  evaluate  the  phase  of  the 
same  function  as  before  for  a  bandwidth  of  f=0.1.  The 
horizontal  line  is  the  log  of  the  magnitude  in  that  region 
with  KR=0.693,  and  the  actual  phase  and  estimated  phase  are 
superimposed  over  each  other.  The  normalizing  constant 
for  the  actual  phase  is  KI=0,104  and  for  the  estimate 


KE=0.082.  The  estimated  phase  is  close  to  the  actual  val¬ 
ues.  Fig. 3*2. 3  shows  the  results  when  the  bandwidth  is 


increased  to  f=0.5*  The  solid  line  is  the  log  of  the 
magnitude  with  KR=0.693»  the  dotted  line  is  the  actual 
phase  with  KI=0.482,  and  the  estimate  of  the  phase  is  given 
by  the  dashed  line  with  KE=0.520.  By  increasing  the  band¬ 
width  of  the  signal  we  see  that  the  estimate  is  less  accu¬ 
rate  because  the  pole  is  having  a  larger  effect  on  the 
results.  In  Fig.  3*2.4  we  have  further  increased  the  band¬ 
width  to  f=2.0  with  the  solid  line  representing  the  log 
magnitude  with  KR=0.693i  the  dotted  line  representing  the 
actual  phase  with  KI=1.125,  and  the  estimate  given  by  the 
dashed  line  with  a  normalizing  constant  KE=0.473.  We  see 
that  here  there  is  a  very  large  error  between  the  estimate 
and  the  actual  value  because  of  the  pole.  In  order  to 
improve  the  results,  methods  must  be  introduced  that  will 
reduce  the  effect  of  this  pole. 

MATHEMATICAL  FILTERING  METHOD 


With  the  method  described  by  Nikajama  and  Asakura  [  46]  we 
will  attempt  to  obtain  a  better  phase  estimate  by  first 
preprocessing  the  magnitude,  using  the  method  where  the 
magnitude  squared  is  smoothed  by  convolving  it  with  a 
filter.  The  filter  we  have  chosen  is  the  Gaussian  filter, 
with  the  understanding  that  there  are  other  filters  which 


could  have  been  used  such  as  the  triangular  or  low-pass 
filter.  The  Gaussian  filter  used  has  zero  mean  and  vari¬ 
ance  c r2.  Its  form  is  given  by 

G(t)=exp(-te  /2  a2  )  (3.2-10) 

In  Fig.  3.2.4  we  saw  that  the  estimate  of  the  phase  is 
very  much  in  error,  and  it  is  in  this  region  that  we  wish 
to  improve  the  estimate.  The  magnitude  was  first  smoothed 
by  the  Gaussian  filter  with  variance  cr2= 4.0.  The  resultant 
log  magnitude  is  shown  in  Fig.  3*2.5  by  a  solid  line  with  a 
normalizing  constant  KR=2.46l.  We  note  that  the  smoothing 
has  raised  the  log  magnitude  above  zero  and  has  made  the 
function  constant  over  a  larger  portion  before  it  began  to 
roll  off  at  the  edges.  The  actual  phase  is  given  by  the 
dotted  line  with  a  normalizing  constant  KI=1.125,  and  the 
estimated  phase  given  by  the  dashed  line  has  KE=2.390.  We 
see  that  the  estimated  phase  is  beginning  to  approach  the 
actual  phase  in  its  shape.  In  our  next  attempt  the  vari¬ 
ance  of  the  filter  is  lowered  to  crz  =  0.5.  The  log  mag¬ 
nitude  after  the  smoothing  is  given  by  the  solid  line  in 
Fig.  3.2.6  with  KR=2.133.  The  actual  phase  is  given  by 
the  dotted  line  having  the  same  KI  as  before.  The  dashed 
line  which  is  the  estimate  of  the  phase  has  a  normalizing 
constant  of  KE=1.981.  The  shape  of  the  estimate  has  not 
changed  from  the  previous  case,  but  the  constant  has  de- 
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creased  towards  the  desired  value.  After  changing  the 
variance  to  cr*  =  o.01  the  estimated  phase  is  very  close  to 
the  actual  phase.  Fig.  3.2.7  shows  the  log  magnitude  giv¬ 
en  by  the  solid  line  with  KR=l.l6l.  The  estimated  phase 
given  by  the  dashed  curve  is  seen  to  follow  closely  the 
shape  of  the  actual  phase  curve  given  by  the  dotted  line. 
The  normalizing  constant  for  the  estimate  is  KE=0.821  as 
compared  to  the  constant  of  the  actual  phase  given  by 
KI= 1.125.  The  constants  are  close  enough  that  we  wanted 
to  see  how  the  real  and  imaginary  parts  of  the  spectrum 
using  the  estimated  phase,  compare  to  the  spectrum  when  the 
actual  phase  is  used.  Fig.  3*2.8  show  the  real  part  of 
the  spectrum.  The  dotted  line  is  the  actual  signal  and  the 
dashed  line  is  the  one  obtained  when  the  estimated  phase  is 
used.  In  both  cases  the  normalizing  constant  is  given  by 
7.389.  The  imaginary  part  of  the  spectrum  is  shown  in  Fig. 
3.2.9.  The  dotted  line  is  the  actual  curve  with  a  nor¬ 
malizing  constant  KI=2.97  and  the  dashed  curve  is  the  es¬ 
timate  with  KE=2.46.  In  Fig.  3*2.10  we  show  the  Fourier 
transform  of  the  estimated  spectrum  to  see  how  closely  it 
relates  to  the  actual  exponential  signal.  As  can  be  seen, 
the  dashed  line  which  represents  the  estimate  is  in  close 
agreement  to  the  actual  exponential  signal.  The  method 
just  described  is  used  before  the  estimate  is  obtained. 
Next  we  would  like  to  discuss  a  method  to  use  after  the 
phase  has  been  estimated. 


CORRECTION  FACTOR  METHOD 

In  a  paper  by  Ford  [56]  a  correction  factor  can  be  found 
if  the  assumption  is  made  that 

1  <#>  (w)-B  |<€  w»Xl  (3.2-11) 

What  this  assumes  is  that  for  |w|  >  £l,  the  phase  can  be 
represented  by  a  constant  with  only  a  small  error  in  the 
assumption.  In  the  case  of  the  exponential  signal,  we 
have  seen  that  the  phase  of  the  spectrum  becomes  constant 
as  it  approaches  f=12.  Therefore,  we  will  use  the  same 
method  as  Ford  to  obtain  a  similar  correction  factor. 

Using  the  modified  Hilbert  transform  we  can  write  the 
phase  <f>  (w)  as  the  sum  of  three  integrals  given  by 

(w)=-l_  f  K(w,y)  <£(y)  dy  -J_  f  K(w,y)  ln|A(y)|dy 

*  J-  co  W-n 


-1  f°°K(w,y)<£  (y)  dy  =  I,+  I«+I,  (3-2-12) 

Trjn. 


where  K(w,y)  is  the  kernel  given  by 
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K(w,y)  =  l-cos(w-y)T 


(3.2-13) 


If  we  make  use  of  (3.2.11)  we  can  rewrite  the  integrals  I, 
and  I 3  as 

I|(w)=-B  f  1-cosz  dz  (3*2-l4a) 

7T  *(w+n.)T  ^ 


(3*2- 14a) 


.  (w-Xl)T 

I,(w)=-B  /•  1-cosz  dz 


■co  z 


(3*2-l4b) 


It  can  be  shown  that  because  of  the  odd  symmetry  of  I(  (w) 
and  Is(w) ,  that  (3.2-14)  can  be  written  in  the  form 


I, (w)+I 


»(w)=-B  f  1 

77 '  J  ~ 


-<w-n.)r 

*  l-cos  z  dz 


(3.2-15) 


<w-il)r  Z 


which  is  the  correction  factor  to  b©  added  to  (3*2-3). 
Equation  (3*2-15)  has  no  closed  form  solution  so  that  the 
integral  has  to  be  evaluated  using  numerical  methods.  In 
Fig.  3*2-11  we  have  evaluated  the  integral  with  B=1  and 
T=1  up  to  f=12.0.  In  Fig.  3*2.12  we  have  added  the  correc¬ 
tion  factor  to  the  phase  estimate  where  the  normalizing 
constant  of  the  estimate  is  KE= 1.355  as  compared  with  the 
actual  KI=1.491.  We  see  that  up  to  f=3*75,  we  have  been 
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able  to  almost  have  an  exact  estimate.  This  tells  us  that 
the  correction  factor  is  good  only  in  the  region  where 
there  is  a  constant  slope.  Once  the  slope  tends  to  zero, 
the  factor  is  of  no  use.  This  is  what  is  observed  in  Fig. 
3.2.12,  that  in  the  region  where  the  slope  was  getting 
smaller,  the  correction  factor  was  unable  to  improve  the 
estimate . 


METHOD  OF  PARTITIONING 


The  method  to  be  described  here  is  based  on  partitioning 
the  spectrum  into  small  intervals  and  estimating  the  phase 
for  that  small  interval.  Once  all  the  estimates  have  been 
obtained  for  each  interval,  they  are  combined  to  form  the 
total  estimate.  If  we  were  to  take  a  look  again  at  Fig. 
3.2.2  we  would  see  that  in  the  region  up  to  f=0.1  we  were 
able  to  obtain  a  close  approximation  of  the  phase.  Based 
on  this  we  want  to  break  the  spectrum  into  intervals  whose 
spacing  is  f=0.1.  What  we  intend  to  do  is  to  shift  the 
log  magnitude  of  that  interval  down  to  the  dc,  and  eval¬ 
uate  the  modified  Hilbert  transform  to  obtain  an  estimate 
of  the  phase,  and  afterwards  shift  the  estimate  back  into 
the  proper  interval.  In  Fig.  3*2.13*  we  have  shifted  the 
log  magnitude  of  the  Fourier  transform  of  the  exponential 
signal  between  f=.l  and  f=.2  to  the  dc  and  because  of  the 
overlap  at  f=0,  we  have  made  it  equal  to  zero.  (We  found 
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that  had  we  added  the  overlap  at  f=0,  we  would  have  ob¬ 
tained  the  same  results.)  By  using  the  modified  Hilbert 
transform  we  obtained  an  estimate  whose  slope  of  the  line 
is  close  to  the  actual  value.  The  slope  of  the  estimate  is 
0.081  as  compared  to  0.102  of  the  phase  in  the  interval 
f=.l  to  f=.2.  In  Fig.  3.2.14  we  have  taken  the  interval 
f=2.0  to  f=2.1  and  shifted  it  to  obtain  the  phase.  The 
slope  of  the  estimated  phase  is  0.0?3  and  the  actual  slope 
is  0.019.  We  see  that  the  estimate  is  increasing  more 
rapidly.  Fig.  3.2. 15  shows  the  results  when  the  interval 
taken  was  between  f=6.5  and  f=6.6.  The  slope  of  the  esti¬ 
mate  is  O.O32  as  compared  to  the  slope  of  the  actual  phase 
of  0.002.  We  see  that  the  slope  of  the  estimate  follows 
the  phase  of  the  actual  phase  for  the  lower  frequency  inter¬ 
vals.  In  Fig.  3*2.16  we  have  taken  the  region  to  f=1.0  and 
divided  it  into  10  intervals,  each  of  which  had  its  log 
magnitude  shifted  so  that  the  phase  could  be  determined. 
Once  the  phase  was  determined,  it  was  shifted  back  to  its 
orginal  interval  and  added  to  to  the  previous  results.  We 
have  the  slope  of  the  estimate  given  by  0.804  which  is  very 
close  to  the  actual  slope  of  the  phase  of  0.808.  Ther¬ 
efore,  we  have  that  in  the  case  of  a  spectrum  that  contains 
a  pole  in  the  upper  half  of  the  z-plane,  this  method  has 
given  good  results  in  estimating  the  phase  when  using  the 
modified  Hilbert  transform. 


«  _  « 


PERIODIC  ANALYTIC  SIGNALS 


In  this  section  we  would  like  to  take  a  look  at  a  special 
class  of  signals  known  as  periodic  analytic  signals.  These 
signals  are  periodic  and  have  the  property  that  the  spec¬ 
trum  contains  components  for  positive  frequencies  only. 
This  is  the  dual  of  causal  time  signals.  We  are  given  a 
time  signal  f(t)  that  can  be  written  in  the  form 

n 

f(t)=  II(l-ai  expj(Xlt-0i  ))  (3«2-l6a) 

n 

=  £oc*  exp  jkilt  (3*2-l6b) 

In  (3.2-l6b)  we  recognize  this  to  be  the  Fourier  series 
representation  of  the  signal  f(t).  The  fundamental  fre¬ 
quency  is  given  by  &  and  we  see  that  the  signal  f(t)  has 
a  bandwidth  given  by  ,  and  that  the  signal  is  repre¬ 

sented  in  terms  of  only  positive  frequencies.  Equation 
(3»2-l6a)  is  a  factored  form  of  (3*2-l6b)  from  which  the 
location  of  the  zeros  can  be  obtained.  The  reason  why  the 
zeros  of  f(t)  are  of  interest  is  because  in  order  to  obtain 
the  phase  of  f(t)  from  its  magnitude,  all  of  the  zeros  must 
lie  in  the  lower  half  of  the  complex  z-plane.  Once  the 
signal  has  been  factored,  all  one  has  to  do  is  to  check  the 
magnitude  of  a,  .  If  the  magnitude  of  a  is  less  than  one, 
then  the  zero  of  that  term,  zi,  is  in  the  lower  half 
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plane.  To  see  just  where  the  zeros  are  located,  we  first 
consider  when  a  is  positive  with  a  magnitude  less  than 
one.  Then  the  zeros  are  located  at 

z  =  2  vn+Q,  +  j  In  |a,|  (3.2-17) 

SI  ~~si 

where  we  see  that  since  the  magnitude  is  less  than  one,  the 
log  will  be  negative  and  the  zeros  are  in  the  lower  half 
plane.  We  notice  that  the  location  of  the  zeros  are  peri¬ 
odic.  If  the  value  of  a  is  negative  and  the  magnitude  is 
less  than  one,  the  zeros  are  located  at 

Z|  =  2ir(n+l)  +  &i  +  j  In  jajj  (3.2-18) 

SI  SI 

again  showing  that  the  zeros  are  periodic  and  located  in 
the  lower  half  of  the  z-plane. 

Let  us  examine  the  case  when  n=l.  We  have  that  the  signal 
f(t)  is  given  by 

f(t)=l  -  a  expj(ftt-e)  (3.2-19) 


from  which  we  obtain  the  magnitude  and  phase  as 


(3.2-20a) 
( 3 • 2-20b ) 


|f(t)|*=  (l+a*)-2a  cos(£lt-@  ) 

<£(t)  =  arctan  -  a  sin(ftt -6  ) 

1  -  a  cos (il  t-0  ) 

Using  the  magnitude  of  (3.2-20)  and  the  Hilbert  rela¬ 
tionship  of  (3*2-4)  we  want  to  obtain  an  estimate  of  the 
phase  given  by  (3*2-20).  The  constants  that  were  used  are 
a=0.2,  T=0.5  and  f=1.0  with  M=256.  The  results  are  shown 

in  Fig.  3*2.17.  the  solid  line  is  the  log  magnitude  with  a 
normalizing  constant  given  by  0.233*  The  actual  phase  is 
shown  by  the  dotted  line  with  a  normalizing  constant  of 
KI=0.201.  The  estimated  phase  is  shown  by  the  dashed  line 
with  KE=0.201.  We  can  see  that  the  estimate  is  a  good  one 
with  the  graph  almost  identical  to  the  actual  graph. 

In  the  next  example  we  wanted  to  investigate  the  case  when 
n=2  and  the  coefficients  a(  are  both  equal  to  0.2.  Ther¬ 
efore,  we  have 

f(t)  =  (1  -0.2  expU&t ))*  (3*2-21) 

The  bandwith  of  this  signal  is  2  SI  .  For  the  choice  of 
f=l  and  T=12,  the  resulting  phase  is  obtained  in  Fig. 
3.2.18  which  was  obtained  by  taking  the  magnitude  of 
(3.2-21)  and  using  (3.2-4).  The  actual  phase  is  given  by 
the  dotted  line  with  KI=0.403  and  the  dashed  line  is  the 
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estimated  phase  with  KE=0.4l6.  We  again  see  that  the  es¬ 
timated  phase  is  close,  and  all  because  the  zeros  of  the 
function  are  in  the  lower  half  plane,  which  makes  it  possi¬ 
ble  to  use  the  logarithmic  Hilbert  transform  for  the  deter¬ 
mination  of  the  phase. 

In  the  previous  cases  we  were  investigating  the  signal 
given  in  (3*2-16).  We  next  would  like  to  take  a  look  at 
the  following  type  of  signal 

f  (t  )=exp(-j  ilt)  -  a  exp(-j@  )  (3.2-22) 
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One  difference  in  this  form  of  the  signal  is  that  rather 
than  f(t)  having  a  periodic  set  of  zeros,  there  is  only  one 
given  by 

z=J_  +  j  In  la |  (3*2-23) 

Si  Si 

if  the  constant  a  is  positive,  and 
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z=  Q-tt  +  j  In  |a  | 

Si  Si 


(3.2-24) 


if  the  constant  a  is  negative.  If  we  look  at  the  mag¬ 
nitude  and  phase  we  have 


f(t)  =  (l+a2  )-2a  cos(  Sit-Q  ) 


(3.2-25a) 


■  V. 


<p(t)  =  -  arctan  sinilt  -  a  sin© 
cos  lit-  a  cos  0 


(3.2-251)) 


The  first  thing  that  we  notice  is  that  the  magnitude  of 
(3.2-25)  is  identical  to  (3*2-20a)  but  that  the  phases  are 
not  the  same.  By  choosing  |a[<l  we  should  be  able  to 
obtain  the  phase  from  the  magnitude  by  using  the  Hilbert 
transform.  But  we  have  already  seen  that  by  using 
(3.2-25a)  the  phase  that  we  obtain  is  given  by  (3*2-20b) 
and  by  (3»2-25b).  To  see  what  is  happening  we  have  to 
look  at  the  Fourier  transform  of  (3*2-19)  and  (3*2-22). 
The  Fourier  transform  of  (3*2-19)  has  two  impulses  that  are 
on  the  positive  side  of  the  spectrum  and  (3.2-22)  has  two 
impulses  that  are  on  the  negative  side  of  the  spectrum. 
Therefore,  we  see  that  we  are  unable  to  use  the 
phase-magnitude  relationship  for  (3.2-25)  because  of  the 
position  of  the  spectrum  it  occupies.  Equation  (3*2-22) 
can  be  rewritten  in  the  form 


f (t)=exp- j  £lt(  1  -  a  expj(&t-0  )) 


(3*2-26) 


where  the  bracketed  term  is  the  same  as  (3. 2-19).  There¬ 
fore,  to  obtain  the  phase  of  (3*2-25b)  we  must  first  modu¬ 
late  the  signal  of  (3*2-22)  so  that  is  in  the  positive  side 
of  the  spectrum.  After  this  the  phase  of  the  signal  can  be 
obtained  using  the  logarithmic  Hilbert  transform.  Once 
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the  phase  is  found,  the  linear  phase  -  ilt  is  added  so  that 
we  obtain  the  phase  of  (3.2-25). 


i.3  Conclusion 


We  have  shown  a  new  form  of  the  Hilbert  transform,  which 
has  the  property  that  it  possesses  no  singularity.  It  has 
been  applied  to  determining  the  phase  from  the  magnitude  of 
the  special  class  of  signals  that  have  no  zeros  in  the 
upper  half  of  the  z-plane.  To  improve  the  estimates  be¬ 
cause  of  poles,  we  have  used  smoothing,  a  correction  fac¬ 
tor,  and  a  new  method  of  partitioning.  One  advantage  of 
the  new  transform  is  that  in  evaluating  the  convolution  we 
do  not  use  the  fast  Fourier  transform.  The  reason  for  the 
advantage  is  that  in  evaluating  the  convolution  using  the 
fast  Fourier  transform,  the  frequency  spread  of  the  func¬ 
tion  may  exceed  the  frequency  domain  defined  by  the  inverse 
of  the  sampling  distance,  and  as  a  consequence  the  desired 
Fourier  transform  of  the  function  is  distorted.  Thus  the 
evaluation  of  the  convolution  yields  a  large  error. 
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NORMALIZED  AMPLITUDE 


Fig.  3.1.4  Estimated  real  part  using  the  modified 
Hilbert  transform  with  normalizing  constant 
KR=0.5,  KI=0 . 362  and  KE=0.494. 
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Fig.  3-1-5  Estimated  imaginary  part  using  the 
modified  Hilbert  transform  where  TBP  is  reduced 
to  3  with  normalizing  constants  KR=0.5,  KI=0.362 
and  KE= O.356. 


FREQUENCY 


Fig.  3.I.6  Estimate  of  the  imaginary  part  of  the 
spectrum  of  an  exponential  function  using  the 
modified  Hilbert  transform,  with  normalizing 
constants  KR=2.0,  KI=1.0  and  KE=1.0. 


Fig.  3.2.1  Estimate  of  phase  from  magnitude 
measurement  of  the  Fourier  transform  of  the  ex¬ 
ponential  signal  with  normalizing  constants 
KR-1.841,  KI=1 .491  and  KE=2.56l. 
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NORMALIZED  UNITS 


Fig.  3.2.?  Estimated  phase  of  the  spectrum  of  an 
exponential  signal  after  Gaussian  filtering  with 
o*=0.01  where  KR=l.l6l,  KI  =  1.125  and  KE=0.821. 
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Fig. '3»2.l4  Phase  estimate  of  the  shifted  inter¬ 
val  from  f=2.0  to  2,1  of  the  spectrum  of  an  ex¬ 
ponential  signal  where  KI=0.019  and  KE=0.073. 
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Fig.  3.1.15  Phase  estimate  of  the  sh 
val  from  f=6.5  to  6.6  of  the  spectrui 
ponential  signal  where  KI=0.002  and  KE 


Fig..  3.2,16  Estimate  of  the  phase  of  the  spec¬ 
trum  of  the  exponential  signal  after  having  par¬ 
titioned  the  interval  into  ten  segments  with 
KI=0. 808  and  KE=0.804. 
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Chapter  4 


SUMMARY  AND  EXTENSIONS 


4.0  Introduction 


A  new  method  of  resolving  two  frequencies  spaced  closely 
together  buried  in  noise  has  been  presented.  The  method 
involves  knowing  certain  information  about  the  signal  a 
priori.  This  information  was  used  to  give  a  set  of  equa¬ 
tions  where  there  are  less  unknowns  than  there  are  equa¬ 
tions.  Since  the  system  of  equations  is  overdetermined, 
there  are  many  possible  solutions.  We  have  selected  the 
solution  vector  which  yields  the  minimal  norm  of  the  resid¬ 
ual  error.  By  so  doing,  we  were  able  to  resolve  the  two 
frequencies  of  the  signal  beyond  the  limit  imposed  by  the 
uncertainty  principle  of  signal  processing. 


A  second  topic  that  was  discussed  is  a  new  form  of  the 
Hilbert  transform.  An  advantage  of  this  new  form  is  that 
it  takes  into  account  for  signals  that  have  been  observed 
for  short  time  durations.  One  use  of  the  Hilbert  transform 
is  in  phase  estimation  from  magnitude  measurements.  This 
can  be  accomplished  if  the  signal  is  causal  and  has  zeros 
that  occur  only  in  the  lower  half  of  the  complex  z-plane. 
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If  zeros  occur  in  the  upper  half  of  the  z-plane,  methods 
have  to  be  introduced  to  reduce  the  effects  of  the  zeros 
on  the  phase  estimation.  A  new  method  for  estimating  the 
phase  was  shown  where  the  spectrum  is  partitioned.  All 
methods  were  demonstrated  via  computer  simulations. 

The  first  chapter  of  this  thesis  contains  material  of  an 
introductory  nature.  In  Chapter  2,  a  new  method  of  resolv¬ 
ing  two  frequencies  spaced  closely  together  when  the 
time-bandwidth  product  is  very  low  was  presented.  The 
method  was  shown  to  give  favorable  results  under  low 
signal-to-noise  ratios.  In  the  process  of  resolving  the 
frequencies,  we  obtained  an  estimate  of  the  time  signal. 
In  Chapter  3,  we  introduced  a  new  form  of  the  Hilbert  trans¬ 
form.  With  the  new  form  of  this  transform  we  showed  how 
the  phase  can  be  estimated  from  the  magnitude  of  a  signal, 
given  that  certain  conditions  are  met.  When  these  condi¬ 
tions  cannot  be  met,  modifications  have  to  be  performed 
before  or  after  the  estimate  is  obtained. 

4.1  Extensions  of  Methods 

There  are  many  extensions  that  can  be  applied  to  the  meth¬ 
ods  discussed  in  this  thesis.  In  the  case  of  spectral 
estimation,  if  it  is  known  that  there  are  three  or  more 
frequencies  involved,  then  it  is  possible  to  obtain  a  set 


of  equations,  and  constraints  on  the  solution  vector  in  an 
identical  manner  as  was  described  in  Chapter  2.  In 
two-dimensional  signal  processing  such  as  is  found  in 
optics,  When  the  two-dimensional  signal  is  seperable,  it 
is  posible  to  use  this  method  on  each  dimension  to  obtain 
better  resolution  for  impulsive  type  spectra.  This  can 
be  further  investigated  to  be  used  for  image  enhancement 
of  electron-microscope  images,  bandwidth  compressed  video 
images,  and  image  enhancement  of  medical  diagnostic  im¬ 
ages  . 

For  the  modified  Hilbert  transform,  the  formulation  can  be 
extended  to  n-dimensions  if  the  signals  are  separable  in 
each  dimension.  If  involved  with  superresolution,  methods 
can  be  developed  that  will  utilize  some  superresolving 
methods  to  obtain  a  superresolving  Hilbert  transform.  This 
could  be  used  as  a  constraint  for  the  estimation  of „ causal 
signals  that  have  reduced  time-bandwidth  products. 


APPENDIX 


Suppose  we  are  given  a  set  of  data  points 

c(i) ,  i=l ,2,3, ... ,m  (1) 

and  a  model  depending  upon  certain  parameters  x 

d(i,x) ,  i=l ,2,3, . . • ,m  ( 2 ) 

where  x  is  a  column  vector  containing  n  elements  and  we 
have  m>  n.  A  way  of  measuring  how  far  the  model  is  from  the 
data  for  any  choice  of  x  is  defined  as  the  residual  f(x), 
and  is  a  function  of  x,  d(i,x),  and  c(i).  Each  of  the 
c(i)  is  a  real  number,  the  d(i,x)  are  real-valued  functions 
of  x  and  f  is  a  real-valued  nonnegative  function.  The 
data-fitting  problem  consists  of  selecting  the  parameters 
x  so  as  to  minimize  f.  These  problems  have  more  structure 
than  do  general  minimization  problems  because  f  is  often  of 
a  form  such  as 

f (x)=  Z  r(c(i)-d(i,x) ) 

f(x)=max  r(c(i)-d(i,x) )  (3) 

for  some,  simple,  nonnegative  function  r(.). 


This  data-fitting  problem  is  linear  if  d(i,x)  has  the  form 


d(i,x)=Dx  (4) 

where  D  is  an  mxn  matrix.  If  all  x  are  feasible  for  con¬ 
sideration  in  finding  the  minimum  of  f,  the  problem  is 
unconstrained.  If  the  x  are  restricted  to  lie  in  some  set 
S,  then  the  problem  is  constrained.  A  common  way  to  spec¬ 
ify  such  feasible  sets  S  is  through  the  use  of  equalities 
and  inequalities 

p(j.x)=q(j)  ,  3*1,2 . i| 

g(k, x)=h(k)  ,  k=l , 2 . lj  (5) 

As  an  example  of  such  problems  one  mdy  have  (4)  where  f 
involves  summation  with  r(.)=(  )  .  This  leads  to  the 

familiar  linear  least- squares  problem 

min  f(x)=  JJ c— Dx jj  =  £  (c(i)-Dx)2.  (6) 

This  problem  becomes  a  linearly  constrained  problem  if, 
for  example,  the  elements  of  x  must  all  be  nonnegative. 

If  one  was  to  restrict  oneself  to  completely  linear  prob¬ 
lems  in  which  r=  j.j  ,  then  this  would  involve  1»  and  l<x> 
data-fitting.  The  least-squares  (la  )  criterion  for  meas- 


uring  the  fit  of  a  model  to  data  is  the  most  reasonable  one 
to  use  if  it  is  known  that  the  data  c(i)  are  approximations 
to  certain  quantities  c’(i) 


c(i)=c ' (i)+e(i)  (7) 

where  the  errors  e(i)  are  all  normally  distributed  with 
zero  mean  and  common  variance.  This  is  the  case  when  the 
c(i)  result  from  careful,  bias-free  observations  of  well 
behaved  systems.  However,  if  this  is  not  the  case  then  two 
other  situations  are  very  important: 

(1)  the  observations  are  not  always  careful  or 

bias-free,  or  the  system  is  not  always  well-behaved  (1,  ),. 

(2)  the  measurements  are  exact,  and  all  errors  are  zero 
(loo)  • 

In  (1)  there  may  be  occasional  errors  e(i)  which  are  quite 
wild.  It  is  desirable  to  ignore  the  corresponding  observa¬ 
tions  and  fit  the  model  only  to  the  good  data.  In  (2),  it 
is  desired  to  get  the  model  evenly  close  to  every  one  of  the 
observations  c(i).  The  1,  norm  can  be  expressed  as 

f(x)=  £  jc(i)-DxJ  =  JJc-Dxjl  (8) 

«  * 

and  the  1  co  norm  can  be  expressed  as 


f(x)=max|  c(i)-Dx|=  f J c— Dx{ j oo 


(9) 


The  former  case  is  of  interest  to  statisticians  (as  robust 
regression)  and  the  latter  case  to  mathematicians  (function 
approximation) . 

To  see  the  difference  between  the  three  norms,  Fig.  (A.l) 
shows  observations  made  of  a  45°  line  over  a  unit  square. 
If  these  observations  are  subject  to  minor  random  fluctua¬ 
tions,  and  if  an  attempt  is  made  to  fit  the  observations 
with  a  straight  line,  the  lco  »  I2  *  lj  norms  result  in 
approximately  the  same  solution  as  seen  in  Fig.  (A.l).  On 
the  other  hand,  if  one  of  the  observations  is  widely  at 
odds  with  the  others,  then  the  1,  line  will  not  be  dis¬ 
turbed,  the  12  line  will  be  displaced  somewhat  towards 
the  wild  point,  and  the  Iqt)  line  will  be  centered  between 
the  wild  point  and  the  remaining  ones  as  shown  in  Fig. 
(A. 2). 

It  is  for  these  reasons  that  the  solution  to  the  problem  of 
section  3*1  will  be  solved  for  using  the  1,  norm.  The 
data  c(i)  is  very  noisy  causing  wild  fluctuations  in  the 
data,  and  it  is  the  1(  norm  that  will  fit  the  model  to 
the  good  data  only. 

The  linear  problem  of  section  3.1  can  be  written  in  the 
form 


minimize  £  (c(i)-Dx) 

i 

subject  to  Px=q 

Gx=h  (10) 

whose  solution  can  be  arrived  at  using  linear  programming. 
Since  the  solution  x  and  the  residual  need  both  be  negative 
as  well  as  positive,  and  because  linear  programming  gives 
only  positive  solutions,  31”  additive  variables  must  be 
introduced.  This  gives  us 

c-Dx+dy=v-u 

x*0 

y-0 

v-0 

u-0  (11) 


so  that  (10)  can  be  written  as 

minimize  £  u+v 

subject  to  Px-Py  =q 

Gx-Gy  =h 

Dx-Dy-u+v=c  (12) 
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which  can  be  solved  by  a  method  such  as  the  simplex  method 


Fig.  A.l  Observations  of  a  45  line  over  a  unit 
square  showing  the  1,,  12  and  lco  curve  fits  of 
the  line. 
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